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ABSTRACT 


An  orbiting  detector  of  infrared  (IR)  energy  may  be  used  to  detect  the  rocket 
plumes  generated  by  ballistic  missiles  during  the  powered  segment  of  their  trajectory.  By 
measuring  angular  directions  of  the  detections  over  several  observations,  the  trajectory 
properties,  launch  location,  and  impact  area  may  be  estimated  using  a  nonlinear  least- 
squares  iteration  procedure.  Observations  from  two  or  more  sensors  may  be  combined  to 
form  stereoscopic  lines-of-sight  (LOS),  increasing  the  accuracy  of  the  estimation 
algorithm.  The  focus  of  this  research  has  been  to  develop  a  computer-model  of  an 
estimation  algorithm,  and  determine  what  parameter,  or  combination  of  parameters  will 
significantly  affect  on  the  error  of  the  tactical  parameter  estimation.  This  model  is  coded 
in  MATLAB,  and  generates  observation  data,  and  produces  an  estimate  for  time,  position, 
and  heading  at  launch,  at  burnout,  and  calculates  an  impact  time  and  position.  The  effects 
of  time  errors,  LOS  measurement  errors,  and  satellite  position  errors  upon  the  estimation 
accuracy  were  then  determined  using  analytical  and  Monte  Carlo  simulation  techniques. 
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1.  INTRODUCTION 


The  TRW-buUt  Defense  Support  Program  (DSP)  satellites  have  been  the 
spacebome  segment  of  NORAD’s  Tactical  Warning  and  Attack  Assessment  system  since 
the  early  1970’s.  Using  infrared  detectors  that  sense  the  heat  from  missile  plumes  against 
the  Earth  background,  these  orbiting  sentries  detect  ballistic  missile  launches.  DSP  also 
detects  nuclear  detonations  in  support  of  Nuclear  Test  Ban  monitoring.  The  DSP  system 
provides  near  real-time  detection  information  in  support  of  DOD’s  tactical  warning  and 
attack  assessment  mission  and  is  supported  by  a  network  of  fixed  and  mobile  ground 
stations  that  process  and  disseminate  information  to  military  commanders  worldwide.  The 
Cold  War  mission  of  the  DSP  system  was  to  detect  massive  intercontinental  ballistic 
missile  (ICBM)  attacks.  United  States’  response  to  such  an  attack  only  required  timely 
and  unambiguous  warning.  Missile  flight  times  were  much  longer  than  the  time  required 
to  launch  a  retaliatory  attack.  Precise  radar  tracks  could  be  established  with  enough  time 
to  mitigate  effects  as  much  as  technology  allowed. 

The  current  combat  environment  demands  much  more  from  launch  detection 
satellites.  The  present  threat  is  from  tactical  ballistic  missiles  (TBM’s),  which  exhibit 
much  cooler  burn  and  shorter  thrust  times,  and  possibly  more  depressed  trajectories  than 
those  exhibited  by  ICBM’s.  TBM’s  can  be  launched  from  almost  anywhere  within  a  large 
geographical  area  of  interest,  \vith  lofted  or  depressed  trajectories.  There  may  be  no  other 
sensors  to  quantify  impact  parameters  in  time  to  employ  anti-ballistic  missile  (ABM) 
weapons  or  alert  potential  victims  within  the  impact  zone. 

For  budgetary  reasons,  the  United  States  is  forced  to  use  existing  launch  detection 
satellites  to  counter  the  TBM  threat  into  the  beginning  of  the  next  century.  [Ref  1]  More 
rapid  extraction  of  more  precise  information  from  these  existing,  technology-limited 
satellites  must  be  accomplished  in  the  interim  of  acquiring  the  next-generation  TBM 
detection  satellite  system.  Anti-ballistic  missile  (ABM)  systems,  such  as  Patriot  or  Aegis, 
may  be  employed  as  ABM  umbrellas  in  a  tactical  area  of  interest,  provided  they  receive 
precise  and  timely  cueing  from  these  detection  systems. 
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This  raises  the  question;  “How  accurate  is  the  information  provided  by  DSP?”  To 
begin  to  answer  this  question,  an  engineering  error  analysis  must  be  performed  to 
determine; 

•What  are  the  errors  present  in  the  detection  system? 

•  Which  errors  have  the  greatest  effect  upon  the  accuracy  of  DSP  output 
information? 

•What  are  the  effects  of  the  errors  on  the  results  (individually  and  collectively)? 

By  modeling  the  algorithm  used  by  the  tactical  warning  system  to  determine  trajectory 
elements  and  predict  impact  zone  location  of  a  detected  TBM,  it  is  then  possible  to 
introduce  the  inherent  errors  separately  into  the  model,  and  then  study  their  effects  upon 
the  results.  The  relative  magnitudes  of  their  effects  upon  the  output  and  their  effects  upon 
the  results  are  then  determined.  This  is  done  both  qualitatively  and  numerically 
(statistically),  and  the  results  are  compared. 
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II.  BACKGROUND 


The  DSP  system  consists  of  one  or  more  satellites  in  geosynchronous  orbit  and 
one  or  more  ground  receiving  stations.  A  DSP  satellite  consists  of  a  bus  (spacecraft) 
segment  and  a  payload  (sensor)  segment.  The  satellite  is  10  meters  long,  7  meters  in 
diameter,  and  weighs  over  2300  kilograms.  [Ref  2] 


Figure  1.  Diagram  of  DSP  Satellite  [From  Ref  2] 

The  spacecraft  segment  receives  and  transmits  all  commands  and  data,  provides 
attitude,  spin-rate,  and  station-keeping  functions,  and  provides  and  distributes  satellite 
power.  Solar  cells  mounted  on  the  spacecraft’s  cylindrical  body  and  on  four  solar  panels 
mounted  opposite  the  sensor  generate  the  electrical  power.  The  sensor  collects  IR  energy, 
provides  onboard  (initial)  data  processing,  and  supplies  precise  orientation  data. 
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The  satellite  is  placed  in  a  circular,  equatorial,  geosynchronous  orbit  (Titan/IUS), 
oriented  so  that  the  telescope  is  pointed  toward  the  earth,  and  spun  along  its  longitudinal 
axis  at  six  rpm.  This  configuration  is  called  a  “yaw  spinner”.  The  axis  of  rotation  of  the 
satellite  is  normal  to  the  earth’s  surface  (nadir-pointing).  The  satellite’s  spin  allows  the 
sensor  to  regularly  scan  the  earth  and  distributes  the  thermal  load  uniformly.  The  spin  also 
allows  one  set  of  attitude  control  components  to  effect  two-axis  earth  pointing  control 
(roll  and  pitch)  by  time-sharing.  A  counter-rotating  momentum  wheel  keeps  the  net  spin 
momentum  near  zero,  thereby  minimizing  the  gyroscopic  effect  due  to  coupling  between 
the  spin  motion  and  the  orbit  rate  with  minimum  fuel  expenditure. 


The  IR  telescope  is  tilted  from  the  spin  axis,  so  that  the  photo-electric  cell  (PEC) 
array  covers  the  radius  of  the  earth.  As  the  satellite  rotates,  the  entire  surface  of  the  earth 
within  the  FOV  is  scanned  by  the  IR  detector,  shown  in  Figures  2  and  3. 
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Detector 

Module 


Figure  3.  Focal  Plane  Array  Scanning  FOV  [From  Ref,  3] 
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Detection  of  IR  sources  is  accomplished  with  the  telescope  and  PEC-array 
portions  of  the  s«isor.  The  signal  dectronics  onboard  the  satellite  partially  discriminate 
the  targets  from  the  background  and  the  ground  station  completes  the  task.  Lxxjation  of 
the  IR  sources  is  derived  from  orientatimi  of  the  IR  telescope  line-of-sight  relative  to  the 
earth’s  surfece.  Star  sensors  augment  the  sensor  data  for  precise  determination  of  sensor 
pointing  direction. 

The  PEC-array  of  the  IR  detector  is  mounted  with  the  nadir  end  at  the  center  of 
rotation  of  the  telescope  (see  Figure  3).  This  array  contains  over  6000  detector  cells  that 
are  sensitive  to  energy  in  the  infrared  wavelengths.  As  the  PEC  array  scans  the  FOV,  a 
cell  passing  across  an  IR  source  will  generate  a  volt^  with  an  amplitude  proportion^  to 
the  signal  intensity.  This  voltage  signal  is  termed  an  IR  return,  and  is  transmitted  to 
ground  processing  stations  after  amplification  and  background  filtering.  A  simplified 
diagram  of  the  data  distribution  network  is  shown  in  Figure  4. 


Figure  4.  Data  Distribution  Diagram  (Ref  4] 

DSP  has  ground  stations  globally  positioned  to  receive,  process,  and  report 
mission  data  to  users.  The  Satellite  Control  Facility  tracks  the  satellite,  transmits 
commands,  receives  mission  and  satellite  status  data,  and  forwards  the  data  to  the  Data 
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Reduction  Caiter  (DRC)  for  further  processing.  The  DRC  processes  the  data,  extracts 
significant  returns  associated  with  missite  launches  and  nuclear  detonations,  and  generates 
event  messages  for  transmission  to  the  users  over  the  Ground  Communications  Network. 
[Ref  2] 
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in.  THE  ALGORITHM 


{Thomas  Jerardi  has  graciously  given  his  permission  to  adapt  and  modify  his 
unpublished pc^r,  “TALON  SHIELD/ALERT  State  Vector  Estimation",  which  is  the 
basis  for  the  majority  of  this  chapter. 5] 

The  detection  of  a  TBM  launch  is  the  starting  point  for  any  ballistic  missile  de¬ 
fense.  DSP  does  this  by  detecting  the  IR  radiation  emitted  by  the  exhaust  plume  of  a 
launching  missile.  With  detections  by  two  or  more  spacecraft  (stereo  observations),  trian¬ 
gulation  of  lines  of  sight  can  be  used  to  more  accurately  estimate  the  boost-phase  trajec¬ 
tory,  which  is  then  used  to  calculate  launch  position,  state  vector  (position  and  velocity)  at 
missile  engine  burnout,  and  impact  position.  Real-time  knowledge  of  the  launch  position 
allows  targeting  of  the  launcher.  Cueing  for  ABM  weapons  systems,  such  as  Patriot  or 
Aegis,  requires  timely  and  accurate  trajectory  information,  which  can  be  propagated  from 
knowledge  of  the  state  vector  at  burnout.  Impact  time  and  position  data  is  extrapolated 
from  the  state  vector  at  burnout,  and  may  be  used  for  warning  personnel  within  the  target 
area.  Therefore,  the  quality  (accuracy)  of  the  trajectory  estimation  process  is  of  para¬ 
mount  importance  to  Tactical  Ballistic  Missile  Defense  (TBMD).  Understanding  the 
algorithms  and  equations  employed  in  the  estimation  process  is  necessary  to  assess  the 
quality  of  the  estimated  parameters. 

The  tactical  parameter  estimation  process  is  composed  of  several  tasks:  initial 
estimate  of  the  tactical  parameters,  nonlinear  least-squares  estimation  (refinement)  of 
tactical  parameters,  burn-out  time  estimation,  state  vector  generation,  and  impact  point 
calculation.  The  tactical  parameters  are. 

•  To  =  Time  of  launch 

•  L  =Loft 

•  <l>o  =  Launch  point  geodetic  latitude 

•  hi-  Launch  point  geodetic  longitude 

•  ho  =  Launch  point  height  above  WGS-84  ellipsoid 

•oto  =  Flight  trajectory  azimuth  (true  heading) 

These  quantities  will  be  explained  in  more  detail  in  the  following  sections. 
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Observational  data  are  used  to  calculate  initial  estimates  of  tactical  parameters  and 
to  choose  the  appropriate  TBM  profile  from  a  database.  The  profile  trajectory  is  then 
used  to  calculate  expected  observations,  which  are  then  compared  to  the  actual  observa¬ 
tions.  The  differences  are  used  to  calculate  changes  to  the  initial  estimates,  and  a  least- 
squares  iteration  is  run  until  the  differences  between  expected  and  actual  observations  are 
sufficiently  small.  The  result  is  an  accurate  determination  of  the  TBM’s  launch  and  trajec¬ 
tory  parameters.  An  estimate  of  burnout  time  is  then  made,  and  the  calculated  trajectory 
is  extrapolated  to  generate  the  state  vector  at  burnout.  This  allows  the  computation  of 
impact  time  and  position  using  simple  ballistic  trajectory  equations.  Each  of  these  tasks  is 
described  in  detail  in  the  following  sections. 

A.  OBSERVATIONAL  DATA 

The  starting  point  is  the  observational  data,  which  are  measurements  of  IR  radiant 
intensities  (IR  returns)  as  a  function  of  time  from  each  of  the  approximately  6000  focal 
plane  elements  in  the  DSP  satellite’s  sensor.  Attitude  information  (star  sensor  measure¬ 
ments,  jet-firing  data,  momentum  wheel  tachometer  data)  is  also  included  in  the  telemetry 
stream  and  used  to  determine  spacecraft  attitude  history.  The  satellite’s  position  is  de¬ 
termined  by  the  Air  Force  Satellite  Control  Network  (AFSCN).  A  global  perspective  of 
the  observation  geometry  is  shown  in  Figure  5. 


Figure  5.  Observation  Geometry  [From  Ref.  5] 
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The  discrimination  of  TBM  IR  returns  from  earth  background  is  a  complex  proc¬ 
ess.  It  will  be  assumed  to  have  been  effectively  done,  and  that  a  table  of  observations 
corresponding  to  a  single,  boosting  TBM  is  available.  An  example  is  shown  in  Table  1. 


Index 

■SI 

S/C 

m 

Intens. 

Azimuth 

(r»d) 

Elevation 

(rad) 

G.H.A. 

(rad) 

Radius 

(km) 

^^Qll 

T 

s/c 

I 

ft 

gha 

5 

R 

1 

129.36 

1 

14.0 

4.061854 

0.115847 

0.174532 

0 

42164.17 

2 

130.30 

2 

23.0 

2.427020 

0.094741 

1.221730 

0 

42164.17 

3 

135.44 

3 

32.0 

2.062181 

0.142310 

1.832595 

0 

42164.17 

4 

139.36 

1 

29.0 

4.062497 

0.115971 

0.174532 

0 

42164.17 

5 

140.30 

2 

29.0 

2.427264 

0.094753 

1.221730 

0 

42164.17 

6 

145.44 

3 

40.0 

2.061969 

0.142405 

1.832595 

0 

42164.17 

7 

149.36 

1 

49.0 

4.063552 

0.116147 

0.174532 

0 

42164.17 

8 

150.30 

2 

60.0 

2.427660 

0.094746 

1.221730 

0 

42164.17 

9 

155.44 

3 

65.0 

2.061752 

0.142493 

1.832595 

0 

42164.17 

Table  1.  Example  Observational  Data 


The  index,  k,  runs  from  1  to  n,  the  total  number  of  observations  (typically  less  than 
20),  and  is  used  as  a  subscript  for  the  remaining  symbols  to  denote  a  particular  observa¬ 
tion.  Time,  Tk,  is  the  time  of  observation  measured  from  midnight  of  the  day  of  the  obser¬ 
vation,  and  runs  between  0  and  86400  seconds.  Spacecraft  Identification,  S/C  ID,  is 
simply  a  notation  used  to  identify  which  satellite  is  making  which  observation.  The  in¬ 
tensity,  Ik,  is  the  radiant  intensity  of  the  IR  return,  and  is  used  to  select  the  type  of  TBM 
being  detected.  Azimuth  angle,  Pk,  is  the  azimuth  of  the  line  of  sight  of  the  IR  return 
measured  clockwise  from  true  south,  and  has  a  range  of  0  to  2%  radians  (0°  to  360°). 
Elevation  angle,  iik,  is  the  elevation  of  the  return  measured  from  nadir,  and  has  values 
between  0  and  0.175  radians  (0°  to  10°).  Satellite  position  is  given  in  spherical  coordi¬ 
nates.  Greenwich  hour  angle,  ghak,  is  the  angle  between  the  satellite’s  nadir  point  and  the 
prime  meridian,  measured  positive  east.  Declination  angle,  5k,  is  the  angle  above  or  below 
the  equator,  measured  positive  north.  (Greenwich  hour  angle  and  declination  have  been 
adopted  instead  of  sub-satellite  point  longitude  and  latitude  to  avoid  confusion  with  TBM 
latitude  and  longitude.)  Radius,  Rk,  is  the  distance  from  the  satellite  to  the  earth’s  center, 
measured  in  kilometers.  Geosynchronous  satellites  typically  stay  within  a  few  degrees  of 
the  equator,  and  have  a  radius  of  approximately  42,164.17  km. 
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B.  TBM  PROFILE 


A  TBM  profile  is  a  description  of  the  nominal  powered  flight  trajectory  of  the 


given  TBM.  An  example  TBM  profile  is  given  in  Table  2. 


ISI 

Intensity 

Altitude 

(km) 

Range 

(km) 

Time 

(sec) 

Intensity 

Altitude 

(km) 

Range 

(km) 

0 

0.000 

0.000 

32 

60.6 

7.023 

3.195 

1 

■esh 

0.006 

0.000 

33 

62.4 

7.469 

3.491 

2 

36.6 

0.026 

0.000 

34 

64.2 

7.928 

3.803 

3 

36.9 

0.058 

0.000 

35 

66.0 

8.402 

4.132 

4 

37.2 

0.103 

0.000 

36 

68.4 

8.890 

4.479 

5 

37.5 

0.163 

0.001 

37 

70.8 

9.393 

4.844 

6 

37.8 

0.235 

0.004 

38 

73.2 

9.911 

5.229 

7 

38.1 

0.322 

0.010 

39 

75.6 

10.444 

5.633 

S 

38.4 

0.423 

0.020 

40 

78.0 

10.992 

6.057 

9 

38.7 

0.537 

0.036 

41 

81.2 

11.556 

6.502 

10 

39.0 

0.666 

0.058 

42 

84.4 

12.136 

6.969 

11 

39.5 

0.809 

0.087 

43 

87.6 

12.732 

7.459 

12 

40.0 

0.965 

0.124 

44 

90.8 

13.345 

7.973 

13 

40.5 

1.136 

0.171 

45 

94.0 

13.975 

8.511 

14 

41.0 

1.321 

0.226 

46 

96.0 

14.622 

9.075 

15 

41.5 

1.520 

0.292 

47 

98.0 

15.288 

9.665 

16 

42.0 

1.733 

0.367 

48 

100.0 

15.972 

10.282 

17 

42.5 

1.962 

0.453 

49 

102.0 

16.675 

10.928 

18 

43.0 

2.204 

0.550 

50 

104.0 

17.397 

11.604 

19 

43.5 

2.460 

0.658 

51 

104.6  1 

18.140 

12.309 

20 

44.0 

2.731 

0.777 

52 

105.2 

18.904 

13.045 

21 

45.0 

3.015 

0.908 

53 

105.8 

19.690 

13.813 

22 

46.0 

3.312 

1.050 

54 

106.4 

20.499 

14.613 

23 

47.0 

3.623 

1.205 

55 

107.0 

21.332 

15.446 

24 

48.0 

3.948 

1.372 

56 

106.4 

22.190 

16.314 

25 

49.0 

4.286 

1.551 

57 

105.8 

23.075 

1.217 

26 

50.6 

4.637 

1.744 

58 

105.2 

23.986 

18.155 

27 

52.2 

5.001 

1.950 

59 

104.6 

24.925 

19.131 

28 

53.8 

5.378 

2.170 

60 

104.0 

25.894 

20.145 

29 

55.4 

5.769 

2.404 

61 

98.0 

26.894 

21.199 

30 

57.0 

6.174 

2.652 

62 

80.0  I 

27.925 

22.293 

31 

58.8 

6.591 

2.916 

62.5 

20.0 

28.450 

22.850 

Table  2.  Sample  TBM  Profile  [Ref  5] 
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A  profile  consists  of  IR  intensity  as  a  function  of  time,  nominal  (maximum  range 
trajectory)  vertical  and  horizontal  ranges  from  the  launch  point  as  functions  of  time,  and 
maximum  bum  time,  W  (62.5  seconds  in  the  example).  The  detected  radiant  intensity 
over  time  is  compared  to  TBM  profiles  in  a  data  base,  and  the  best  match  is  selected  as  the 
type  of  TBM  being  observed.  This  selection  process,  called  “typing”,  is  complex,  and  is 
also  assumed  to  have  been  done  correctly.  A  more  computationally  efficient  (but  equiva¬ 
lent)  representation  of  the  nominal  downrange  distance  and  altitude  profiles  can  be  found 
by  fitting  quartic  polynomials  to  the  data  in  the  profile  tables.  Use  of  the  polynomial  repre¬ 
sentation  during  calculations  precludes  table  look-up  and  interpolation  problems  for  non¬ 
integer  times-of-flight.  The  profile  distance  polynomial  has  the  fonn: 

dp  =ao +a,t  +  a2t^ +a3t^ +a4t'*  (1) 

and  the  profile  height  polynomial  is; 

hp  =bo +b,t  +  b2t^ +b3t^ +b4t‘'  (2) 

where  t  is  time  of  flight,  and  ai..4  and  bi..4  are  coefficients  of  the  fourth-order  polynomial. 
These  coefficients  are  not  computed  within  this  algorithm,  but  are  assumed  to  be  known, 
exact,  correct,  and  available  fi’om  the  typing  process;  i.e.,  perfect  a  priori  knowledge  of 
the  particular  observed  TBM’s  nominal  trajectory.  Since  real-life  TBM’s  do  not  fly  the 
profile  exactly,  using  the  profile  to  determine  the  trajectory  introduces  an  error  into  the 
algorithm.  The  inherent  error  of  inexact  profile  information  is  ignored  in  this  analysis. 
Assuming  the  profile  to  be  exact  does  not  invalidate  the  error  analysis,  however. 

A  “loft”  parameter,  L,  is  used  to  account  for  trajectories  above  or  below  the  nomi¬ 
nal  profile  This  \&y  simple  model  for  a  ballistic  missile  trajectory  adjusts  the  profile  to 
give  actual  range; 

d  =  (l-1.5L)dp  (3) 

and  altitude; 

h  =  (l  +  L)hp  (4) 

Loft  varies  over  approximately  ±0.25.  Figure  6  is  a  plot  of  nominal  (L  =  0),  lofted  (L  = 
+0.25)  and  depressed  (L  =  -0.25)  trajectories. 
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Figure  6.  Tr^ectory  Profiles  for  Three  Loft  Values 

C.  LINE  OF  SIGHT  PROJECTION 

The  least-squares  iteration  process  begins  with  an  initial  estimate  of  the  tactical 
parameters.  In  order  to  make  an  estimate  for  the  initial  TBM  position  (Ao,  <|)o),  the  focal 
plane  measurements  must  be  changed  into  a  LOS  between  the  satellites  and  the  TBM, 
which  points  to  the  position  of  the  TBM  on  the  globe,  The  first  step  is  to  transform  the 
satellite  positions  (gha,  8,  R)  into  earth-centered  rotating  coordinates  (XYZ),  in  which  the 
X-axis  extends  through  the  prime  meridian  and  the  Z  axis  is  aligned  with  earth’s  spin  axis, 
depicted  in  Figure  7. 
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Figure  7.  Coordinate  Frame  Illustration 
The  satellite  position  transformation  from  (gha,  6,  R)  to  (XYZ)  is; 


cos(6^ )  cos(gha^  f 

ys. 

Rjt  cos(6^)sin(ghak) 

UsJ 

1  RkSin(5k)  y 

where,  the  subscript,  s,  denotes  satellite  position,  and  the  sub-subscript  refers  to  the  k* 
observation.  The  focal  plane  coordinates,  r|k  and  Pk,  are  first  transformed  to  the  satellite’s 
Up-East-North  reference  frame,  (UEN),  and  then  must  also  be  transformed  into  (XYZ). 
The  focal  plane  coordinates,  can  be  visualized  with  Figure  8. 
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The  transformation  of  polar  (rik,  Pk)  coordinates  into  Cartesian  (UEN)  coordinates 


is: 


-cos(iik) 

e.,  =  -sin(pj)sin(Tik)  (6) 

-cos(Pk)sin(Tik). 

This  is  the  unit  line-of-sight  (LOS)  direction  vector  in  (UEN)  coordinates.  Next  transform 
into  (XYZ)  (refer  to  Figure  7): 


(7) 


For  the  initial  position  estimate,  it  is  sufficient  to  assume  that  the  line  of  sight 
terminates  on  the  surface  of  a  spherical  earth  with  an  effective  radius  of  reff  =  6371  km. 
The  Cartesian  coordinates  of  the  TBM  are: 


'cos(ghak) 

-sin(ghak) 

o' 

cos(8k) 

0 

-sin(8,,) 

K  = 

sin(ghak) 

cos(ghaJ 

0 

0 

1 

0 

0 

0 

1 

sin(8k) 

0 

cos(8j ) 

~  Pk®k 


Z 


Tk 


(8) 


where  the  subscript,  t,  denotes  TBM  position,  sub-subscript  k  refers  to  the  k'^'  observa¬ 
tion,  and  p  is  the  length  of  the  LOS  vector  as  shown  in  Figure  7.  The  Law  of  Sines  is  one 
of  several  methods  that  may  be  used  to  calculate  p.  The  geometry  of  the  problem  is 
shown  in  Figure  9. 


DSP 


Figure  9.  Initial  Position  Estimate  Geometry 
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The  angle  opposite  p  is: 


r\ 


-1 


RjsinOik)' 


Xk  =7C-'nk-sin 

V  ; 

With  some  trigonometric  identities  and  algebra,  the  expression  for  p  becomes; 

-if|Rk|sin(Tik)^ 


Pk 


_reffSin(Xk)_ 


r.ff  sm 


+  sin 


sin(Tik) 
The  longitude  of  this  point  is; 


sin(Ti^) 


A,k  =  tan 


-1 


and  the  latitude  is; 


(9) 


(10) 


(11) 


(t>k  =  tan 


-1 


(12) 


D.  INITIAL  ESTIMATES 

The  core  of  the  estimation  process  is  a  nonlinear  least-squares  iterative  calculation 
of  the  tactical  parameters.  This  method  is  based  on  a  one-term  (linear)  Taylor  expansion 
of  the  focal  plane  observations  in  terms  of  the  six  tactical  parameters.  This  procedure 
requires  an  initial  estimate  of  the  tactical  parameters.  The  accuracy  of  the  initial  estimate 
does  not  effect  the  accuracy  of  the  final  estimate,  but  can  effect  the  convergence  time  of 
the  least-squares  process. 

The  initial  estimate  of  the  launch  position  and  heading  is  based  on  the  projected 
positions  (cjik  and  X*)  obtained  from  the  LOS  observations.  The  lines-of-sight  from  each 
satellite  usually  do  not  intersect  at  a  single  point  (even  if  they  were  simultaneous  observa¬ 
tions),  as  shown  in  Figure  10.  The  initial  latitude  and  longitude  guesses  are  found  by 
calculating  the  average  position  from  the  first  LOS  observation  from  each  satellite.  The 
assumption  is  made  that  if  one  satellite  can  see  the  TBM,  then  all  can  observe  it.  This 
artificiality  is  not  true  in  the  physical  world. 
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Figure  10.  LOS  Projection  Scatter 

In  this  example,  three  satellites  observe  the  TBM,  so  the  first  LOS  projected  lati¬ 
tudes  and  longitudes  (from  each  satellite)  are  averaged: 

,  (0)  _  (|>|  ~Kt>3 


(0)  _  +  A,2  +X3 

5 - 


(13) 


(14) 


The  superscript,  denotes  the  initial  estimate.  The  last  LOS  projected  positions  are  also 
averaged  to  obtain  the  approximate  final  position  observed: 


Usl  ^ 


(15) 

(16) 
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Equations  (13)  to  (16)  assume  that  the  first  and  last  three  observations  are  from  three 
diflerent  satellites.  If  only  two  satellites  observe  the  launch,  only  the  first  and  last  two 
observations  would  be  averaged. 

The  direction  of  the  last  position  fi'om  the  first  is  the  initial  estimate  for  the  head¬ 
ing, 

a.'"’  -C)cos((l,  -C’l  (17) 

The  function  tan2'’  is  the  two-argument  arctangent,  used  here  because  its  range  is  ±it. 
This  function  is  used  because  a  four-quadrant  answer  is  required. 

The  first  estimate  for  To  is  twenty  seconds  before  the  first  observation: 

T^"^  =  T,  -20  (18) 

and  the  initial  guess  for  L  and  ho  are  both  zero: 

=0  (19) 

h„'"^  =  0  (20) 

E.  EXPECTED  POSITIONS  ALONG  THE  TRAJECTORY 

Given  the  sbc  (estimated)  tactical  parameters  at  launch,  the  expected  position  of 
the  TBM  along  its  boost  trajectory  can  be  determined  by  applying  the  TBM  profile  to  any 


observation  time,  and  in  particular,  Tk,  by  first  computing; 

t,=T,-T,  (21) 

Equations  (1)  through  (4)  then  give  the  expected  altitude  (hk)  and  range  (dk)  fi’om  the 
estimated  launch  point  when  t  =  tk: 

dp|^  ■“  '*’®2fk  "^^sfk  '*'^4fk  (1) 

hp.  =b.+b,n+b,u“+bA'+bA‘  (2) 

dp=(l.|.5L)dp,  (3) 

h,  =(1  +  L)h,_  (4) 
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Next,  it  is  necessary  to  calculate  the  “earth-central  angle”  a  quantity  used  in  (23)  and 
(24): 

=  —  (22) 

^eff 

Ok  is  simply  the  angle  between  the  launch  point  and  the  position  of  the  TBM  with  the 
earth’s  center  as  theyertex.  0  can  be  visualized  with  Figure  1 1 . 


Launch 

position 


Figure  1 1 .  Earth-central  Angle 

Using  0k,  spherical  trigonometry  gives  the  equations  for  the  geodetic  latitude; 


‘l>k  =  --cos'[cos(0Jsin((|>p)  +  sin(0Jcos(<t»o)cos(ao)] 


and  longitude: 


“  ^0  +sin' 


sin(0^)sin(aj) 


(23) 


(24) 


COS((|)J 

The  altitude  is  simply  the  present  height  of  the  TBM  added  to  the  initial  height  at  launch; 

altt=ho+hj  (25) 

In  addition  to  the  geodetic  coordinates,  the  Cartesian  coordinates  are  also  required.  Using 

1 


re  =  6378.137  km,  and  f  = 


298.257 

the  geocentric  latitude  of  the  TBM  is  calculated; 

f  ^  =  tan-’[(l  -  f)^  tan((|)  Jl 


,  the  WGS-84  values  for  the  flattening  of  the  earth. 


(26) 
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as  well  as  the  local  radius: 


Tbo.!,  =  I  (27) 

V(1  -  f)^  cos^(<l>'k  )  +  sin^((|)'^  ) 

These  values  are  then  transformed  into  earth-centered,  Cartesian  coordinates  of  the  TBM; 

cos((|)\  )  +  altk  cosCct)^  )]cos(^^ )  (28) 

Yt,  =  [rfeci,  cos((|)'k  )  +  alt k  cosCit)^ )]  sin(Xt )  (29) 

Zt,  =  fioci,  sin(<|)V  )  +  alt^  sm(^^ )  (30) 


F.  NONLINEAR  LEAST  SQUARES  ESTIMATION 

Since  the  (polar)  focal  plane  coordinates  of  the  observations  are  nonlinear  func¬ 
tions  of  the  tactical  parameters,  a  nonlinear  least-squares  process  is  required.  This  is 
achieved  by  a  sequence  of  linear  least-squares  estimations.  Each  of  these  linear  least 
squares  estimations  is  based  on  a  one-term  (linear)  Taylor  series  expansion  of  the  obser¬ 
vations  in  terms  of  the  tactical  parameters.  A  general  requirement  of  ordinary  linear  least 
squares  is  that  the  noise  contaminating  the  observations  should  be  independent  and  have  a 
common  (preferably  Gaussian)  distribution. 

In  this  case,  the  “polar”  focal  plane  observations  (Pk,  %)  are  transformed  into 
“Cartesian-like”  coordinates  (uk,  Vk); 

Uk  =-tan(TiJsin(p^)  (31) 

v^  =-tan(Ti,)cos(P,)  (32) 

so  that  the  coordinates  have  similar  behavior  with  respect  to  errors  and  noise.  Either  q  or 
sin(q)  could  be  used  in  place  of  tan(q),  but  since  q<10°,  all  three  behave  similarly.  The 
chosen  form  corresponds  to  a  gnomonic  projection  of  the  globe  onto  a  plane  tangent  at 
the  nadir  point.  Figure  12  illustrates  this. 
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Figure  12.  Focal  Plane  Transformation 

Next,  it  is  necessary  to  determine  what  the  focal  plane  coordinates  would  be  for 
each  observation  if  the  TBM  were  actually  located  at  the  predicted  coordinates  given  in 
(28)  to  (30). 

The  expected  line-of-sight  is  computed; 


fAx," 

f  \ 

Ay, 

— 

yr,  -ys, 

IazJ 

-  Zj  j 

In  order  to  determine  the  focal  plane  coordinates  of  azimuth,  P,  and  elevation,  r|, 
the  LOS  vector  given  in  (33)  must  be  rotated  into  the  observing  satellite’s  local  vertical 
coordinate  frame,  (UEN): 

cos(6i)  0  sin(6^)Y  cos(gha^)  sin(ghat)  OVAx^^ 


fu,] 

f 

E, 

InJ 

1 

0  10  I  -sin(ghaj)  cos(gha^)  0 

_  I  0  0  lAAzJ 

This  is  the  transpose  of  the  transformation  matrices  given  in  (7). 

Now  u  and  v  can  be  simply  expressed  in  (UEN)  terms  by  noting  that; 

-E, 


Ay, 


sin(pj  = 
cos(Pk)  = 


Ve,’  +n.’ 

-Nt 


(34) 


(35) 

(36) 
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Expected  focal  plane  coordinates,  (u^ ,  ) ,  are  then: 


(37) 

(38) 


Expected  focal  plane  coordinates  are  denoted  with  a  “  ^  ”,  Taking  the  difference  between 
actual  observations  and  expected  observations  generates  5u  and  5v  (found  on  the  left-hand 
sides  of  (41)  and  (42)): 


8Uk=u^-Uj  (39) 

K=Vk-Vk  (40) 

Now  that  the  changes  in  the  focal  plane  coordinates  are  known,  they  can  be  related  to 
changes  in  the  tactical  parameters.  The  one-term  Taylor  expansions  of  u  and  v  in  terms  of 
the  tactical  parameters  are: 


s:..  _  ^  XT’  ,  ^  ST  .  ^  Cl  ^  ^  Cl  ^  c 

Su  -  ^^0  +  +  — 8<t)o  +  — 5ko  +  — 8h(,  +— 6a„ 

S,._  ^  ST'  ,  ^  SI  .  ^  S-^  ^  S,  ^  c 

ov  -  — 6To  +  — 6L  +  -— 6<|)o  +  — 6ho  +-— Sa,, 

oTq  dL  5(t»o  dXo  5hp 


(41) 

(42) 


These  equations,  written  in  vector-matrix  format,  are  the  linear  least-squares 
model  that  forms  the  basis  of  the  tactical  parameter  estimation  process: 


( 

aij 

5u^ 

au. 

auk"! 

^5To> 

6L 

OT. 

dL 

d^o 

axo 

ah„ 

aap 

5(t>o 

Sv, 

avk 

av. 

avk 

5?io 

1^0 

dL 

^0 

axo 

aho 

aaoJ 

5ho 

¥ 

_ 

(43) 


There  are  two  rows  of  the  center  matrix  in  (43)  for  each  of  the  k  =  (1  to  n)  observations 
for  a  total  of  2n  rows.  This  matrix  of  partials  is  denoted  the  “A”  matrix,  and  represents 
changes  in  the  focal  plane  coordinates  with  respect  to  changes  in  tactical  parameters. 
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The  differences  in  the  focal  plane  coordinates,  5u  and  5v,  are  used  to  generate 
adjustments  to  the  tactical  parameters  by  solving  (43)  for  the  tactical  parameter  variations; 


6L 


Sh„ 

l&XoJ 


The  middle  term  in  (44)  is  the  pseudoinverse  of  A,  and  is  used  because  A  is  an  (2n  x  6) 
matrix,  and  is  generally  not  square.  The  changes  to  the  tactical  parameters,  the  left-hand 
matrix  in  (44),  are  added  to  the  initial  estimate,  and  the  process  is  repeated  until  all  the 
components  of  the  left-hand  vector  in  (44)  are  sufficiently  small,  or  the  convergence  crite¬ 
rion  is  met.  The  result  is  the  correct  tactical  parameters  for  the  observed  TBM  launch. 

The  matrix  of  partial  derivatives.  A,  is  the  key  to  the  whole  iteration  process.  The 
development  of  this  matrix  involves  multiple  use  of  the  Chain  Rule  for  derivatives.  The 
partials  in  (41)  can  be  expanded  to  yield: 


5u_(5u5(|)  dis  dk  du  ^ 

5To  0(|>  UIq  dk  0To  5h  cffo 

du  dk  5u5h 
5L  dL  dkdL  db.dL 

^  ^  3(|>  da  dk  da  dh 

0(|>O  04)  04>o  ^  ^0  ^  ^0 

^  0U  0i|>  0u  0X  ^  0h 

-  — - 1 - 1 - 

dkf,  04*  ^0  ^  ^  ^0 

0u  0u  0J)  da  dk  di 

- = - —  + - + - 

0hp  04)  dk  0ho  da  0ho 

0u  0u  04)  0U  0^  0u  0h 

- = - —  + - + - 

0ap  04)  0ap  dkdUg  0h  0ap 


(45) 

(46) 

(47) 

(48) 

(49) 

(50) 
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In  a  similar  manner,  the  partials  of  (42)  are  expanded; 


dw 

dv  d^ 

dv  dk 

dv  ah 

1 

a^ar. 

1 

ab  ar. 

dv  cl(|> 

dv  &K  av  ah 

Sit)  8L 

ax,  oL  ah  aL 

dw 

_  dv  d^ 

av  ax. 

dv  ah 

1 

ax,  a(j>, 

aha(t>. 

dv 

_  av  5i|> 

dv  ax. 

dv  ah 

dK 

"  okj)  a^o 

ax,  ax,. 

1 

ah  ax. 

dv 

_av  at» 

av  ax. 

dv  ah 

ah„ 

a(i)  ah, 

1 

ax,  ah. 

J 

dti  ah. 

dv 

dv  a(|) 

av  ax. 

av  ah 

1 

da^ 

a(|)  aa. 

j 

ax,  aa, 

,  ah  aa. 

(51) 

(52) 

(53) 

(54) 

(55) 

(56) 


Patterns  and  structure  in  these  twelve  equations  are  more  easily  recognized  when 
these  equations  are  written  in  matrix  form: 


r  au^ 

■  d^ 

ax. 

ah  ■ 

ar. 

ar. 

ar. 

au 

d^ 

ax. 

ah 

dL 

dL 

aL 

(—] 

du 

at> 

dk 

ah 

d^ 

a(|)o 

^0 

d^o 

du 

du 

dk 

ah 

'K 

K 

a?i, 

dk. 

dk. 

ail 

du 

d^ 

ax. 

ah 

ah. 

ah, 

ah. 

ah. 

du 

a|) 

dk 

ah 

\daj 

aa. 

aa. 

aa„_ 
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9^)  5X,  ^ 

ai,  aT„  m, 

^  a  ^ 

SL  cL  5L 

dK  dh 

^0  ^0  ^0 

d^  dk  dh 

dk^  dk^  dko 

d^  dk  db 

aio  cho  chg 

at>  a  ah 

aa^  aa^  aa„ 

Note  here  that  (57)  and  (58)  have  identical  structure.  The  left-hand  sides  are  the  deriva¬ 
tives  of  the  focal  plane  coordinates  with  respect  to  the  tactical  parameters,  the  vectors  on 
the  right  are  the  derivatives  of  the  corresponding  local  plane  coordinates  with  respect  to 
IBM  position  (latitude,  longitude,  and  height),  and  the  center  matrix  describes  the  deriva¬ 
tives  of  the  position  with  respect  the  the  tactical  parameters.  This  effectively  separates  the 
change  in  the  focal  plane  coordinates  with  respect  to  tactical  parameters  into  two  parts:  a 
matrix  that  describes  changes  in  position  due  to  trajectory,  and  a  vector  that  describes  the 
relation  betweai  focal  plane  observations  and  TBM  position. 

The  matrix  is  given  a  name,  Ai: 


5(1) 

dk 

dh 

oTo 

5T, 

dT, 

5(1) 

dk 

dh 

5L 

5L 

ai> 

dk 

5h 

d(f>0 

^0 

d(j> 

dk 

5h 

dk„ 

dk. 

5?.o 

54) 

dk 

5h 

ah. 

dh. 

54) 

dk 

dh 

5a, 

da. 

da, 

(59) 
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so  that  (57)  and  (58)  can  be  rewritten; 


I  aho 
'\daoJ 


The  elements  of  A,  can  be  found  qualitatively  This  is  done  in  Appendix  A,  and 


the  results  are. 


0d  cos((Xo)  0d  sin(cx p ) 

cTo  ~  cT,  cos(4>o) 

13dpcos(ao)  1.5dp  sin(ae) 

rrfcos(4>o) 

1  0 

0  1 

0  0 

dsin(ao)  dcos(ao) 

•"eff  reffCOS(4>o) 


ah 

aio 

hp 
0 

0  I 
1 
0 


(62) 
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The  values  of  the  elements  in  (62)  are  approximations,  but  ^proximate  magni¬ 
tudes  and  correct  signs  are  all  that  are  required  for  the  iteration  process  to  converge.  The 
benefit  of  making  this  simplification  for  A|  is  the  reduction  in  computation  time.  The 
complete  and  exact  formulae  can  be  extrtqjolated  fi’om  Appendix  A. 

Similarly,  vectors  on  the  right-hand  sides  of  (57)  and  (58)  can  be  separated  into 
simpler  components; 


^  _  0u  5x  dady  d[i  dz 
0(j)  dK  dy  dzd^ 

5u  _  5x  ^  dz 

d\  ck  dy  &k  dz  dk 

du  _d\idK  du  dy  du  dz 

5h  ^  0h  0y  ^  Sz  c5h 

dv  _  dv  dK  5v?y  dv  dz 

d^  t?x  3i|)  ^  dz  d^ 

dw  _  dv  d^  d\  dy  dw  dz 

ok  d>i  ck  dy  dk  dz  dk 

dw  _ds/  dK  d\  dy  d\  dz 
dh  d)i  (h  dz  ^h 


(63) 

(64) 

(65) 

(66) 

(67) 

(68) 


As  before,  these  equations  can  be  written  in  matrix  form: 


dx 

dy 

dz 

(du'' 

at) 

0U 

dx 

dy 

dz 

au 

dk 

d?. 

dk 

a 

dy 

cri 

dy 

dz 

au 

c3h 

a 

vaz^ 

^  dv^ 

dx 

dy 

dz 

^av^ 

d^ 

at) 

dx 

5v 

dx 

ay 

az 

dv 

dk 

a 

dk 

a 

dy 

dv 

ck 

dy 

dz 

dv 

a 

ah 

a 

(69) 


(70) 
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The  center  matrix  is  renamed  Az. 


Ox 

dy 

Oz 

0«t) 

Oit> 

Ox 

Oy 

Oz 

a 

dX 

a 

dx. 

dy 

dz 

Oh 

Oh 

Oh 

(71) 


and  its  elements  can  be  approximated  by; 


A, 


'-(r^  +  h)sin(4>)cos(?i)  -(r^  +  h)sin(<i»)sin(X) 
-(teff  +  h)cos(<l))sin(?i)  (r^  +  h)cos(<|))cos(>.) 
cos(<l))  cos(>.)  cos(4>)  sin(X.) 


(r^+h)cos(<|>) 

0 

sin((l») 


The  vectors  on  the  right-hand  side  of  (69)  and  (70)  can  be  further  broken  down: 


“ 

onI 

'Ou> 

5x 

Ox 

Ox 

Ox 

OU 

5u 

OU 

OE 

ON 

Ou 

dy 

Oy 

Oy 

Oy 

ON 

Ou 

Oz 

Oz 

Oz  _ 

loN> 

on] 

rov^ 

5x 

Ox 

Ox 

Ox 

ou 

5v 

ON 

Ov 

dy 

Oy 

Oy 

Oy 

OE 

OU 

ON 

Ov 

^dz 

/ 

-  Oz 

Oz 

Oz . 

The  center  matrices  in  (73)  and  (74)  are  idemical,  and  renamed  A,: 


ON. 

Ox 

Ox 

Ox 

ON 

Oy 

Oy 

Oy 

ON 

Oz 

Oz 

Oz 

As,  th  transformation  from  (XYZ)  to  (UEN)  has  been  derived  previously  in  (7): 


A, 


cosCgha^)  -sin(ghaj 
sinCghaJ  cos(ghaJ 
0  0 


0Tcos(5J  0  -sin(6J 
0  0  10 
1  sin(6i)  0  cos(6J 


(72) 


(73) 


(74) 


(75) 


(76) 
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The  final  step  is  to  determine  the  elements  in  the  vectors  on  the  right-hand  sides  of 
(73)  and  (74).  These  are  easily  computed  from  (37)  and  (38); 


5u  _  E 
_  1 

.^  =  0 

aN 


av  _  N 

aE 


(77) 

(78) 

(79) 

(80) 
(81) 


dw 

Now  (57)  and  (58)  can  be  rewritten: 


2 

U 


ar, 


0 

j  au 

I  ^ 

I  au 

I  ^0 

j  au 
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ah, 

au 

yda^J 


=  A  A  A 


u 

0 


(82) 


(83) 
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( 


5To 

'K 

dv 

dw 

5v 


I 


—  A1A2A3 


0 


V  W 


\.daoJ 


(84) 


Equation  (83)  is  used  to  compute  the  odd  rows  (k  =  odd)  of  A  m  (43),  and  (84)  is 

used  to  compute  the  even  rows  (k  =  even)  of  A  in  (43), 

Focal  plane  coordinates,  (Pi,  in),  are  transformed  into  Cartesian-hke  focal  plane 
coordinates,  (ui,  vi)  ax  (Circled  numbers  refer  to  steps  illustrated  in  Figure  13  )  Satellite 
positions  and  initial  estimates  of  tactical  parameters  are  used  to  generate  expected 
(theoretical)  focal  plane  coordinates,  (ii^  ,Vi)®<s-:,  and  the  A  matrix  ®,  (which  is  denoted 

.  _  in  Figure  13)  and  its  pseudoinverse  (§> .  The  differences  of  the 

^  a{To,L,(|>„,Xosho,ao) 

focal  plane  coordinates,  (5uk,  8vk),  (observed  minus  theoretical)  (§)  are  used  to  generate 
adjustments  to  the  tactical  parameter  estimates  by  computing  . 


/ 

'8u,^ 

8v, 

8L 

8U2 

6<t>0  ' 

<.0  1 

=  [A^Ar’A" 

8V2 

8ho 

l8aoy 

^5v„> 

If  all  the  elements  of  the  left-hand  vector  (adjustments  to  the  tactical  parameters)  are  not 
sufficiently  small,  the  current  values  (j-1)  of  the  tactical  parameters  are  updated  ®  to  the 
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(86) 


<  T  \ 

I 

'5To^ 

L 

L 

5L 

•1>0 

Ho 

''^0  1 

8Xo 

ho 

ho 

5ho 

Uo> 

Uo> 

ISaoJ 

If  the  iteration  process  has  converged,  the  iteration  process  is  terminated  and  the  current 
values  of  the  tactical  parameters  are  the  best  estimates  that  this  process  can  generate  This 
entire  process  may  be  visualized  by  a  flowchart  diagram  drawn  in  Figure  1 3 


Figure  13.  Flowchart  Diagram  of  Iteration  Process 


G.  BURNOUT  TIME  ESTIMATION 

The  estimation  of  burnout  time  is  based  on  the  TBM  profile  maximum  bum  time 
(w),  the  last  observation  time  (T,„),  and  the  nerd  potential  obsen-ation  time  (T.„),  had  it 
occurred.  The  maximum  bum  time  according  to  the  profile  is: 


Two  cases  can  occur: 
(a)ifT„^^>T^then 


(88) 
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(b)  if  Tmax  Tncxt  th©!! 


T^-T,„+|(T„-Tu.)  W 

A  figure  and  example  are  given  below  as  a  demonstration: 


Figure  14.  Burnout  Time  Estimation 

For  the  sample  data  given  in  Table  1,  using  the  profile  given  in  Table  2,  tmax  —  62.5 
seconds,  Ti,*,  =  155.44,  and  the  next  potential  observation  is  T„ext  =  159.36.  With  To  = 
109.36,  the  maximum  bum  time  is  Tnnii  =  171.86,  so  case  (a)  applies  and  burnout  time  is 
estimated  at  Tbo  =  157.40. 

H.  STATE  VECTOR  GENERATION 

The  state  vector  completely  defines  the  TBM’s  position  and  velocity,  and  has  six 
elements.  With  the  tactical  parameters  and  burnout  time  estimates,  generating  the  state 
vector  at  burnout  is  done  by  evaluating  the  position  and  velocity  equations  at  the  burnout 
time: 

t  =Tv  -T„  (90) 

*'bo  *bo  *^0  ^  ^ 
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Substitute  t^,  for  t  in  (1)  and  (2),  and  use  (3),  (4),  and  (22)  through  (25)  to  get  K,  h>o. 


and  altbo: 


(1) 

hp_  =  b,  +  b,t,p  +  bjtj.’  +  bjt  J 

(2) 

dpp=(l-I.5L)d,_ 

(3) 

h,„=(l  +  L)hp^ 

(4) 

Obo  -  — 

*eff 

(22) 

4>b.  =  --cos''[cos(eb,)sin(4>o)  +  sin(9fc.)cos(<^o)cos(ao)] 

(23) 

.  ^  rsin(ei,,)sin(ao) 

coa<*pJ  J 

(24) 

a'tbo  =  ho^hk„ 

(25) 

The  bumout  velocity  is  expressed  in  terms  of  speed  (Vbo),  flight  path  angle  (ybo). 

and  heading  (at*): 

v;„  =  >/d^+h^ 

(91) 

I'- =“■’(!] 

(92) 

/cos((}io)sin(ao)^ 

^bo  ~  1  /A  \  1 

V  cos(({>b„)  J 

(93) 

where  d  and  b  are  the  time  derivatives  of  d  and  h; 

d=^ 

dt 

(94) 

a 

(95) 
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The  state  vector  at  burnout,  Tbo,  is: 

hbo  (96) 

Ybo 

VaJ 

L  CALCULATION  OF  IMPACT  POSITION  AND  TIME 

j  The  equattom  and  methods  in  Ms  section  are  adapted  directly  from  Chapter  6 
of  Bate  et  al  ]  [Ref  6]  The  balSstic  trajectory  is  modeled  in  three  phases;  powered 
flight,  which  is  the  portion  tha,  DSP  observes,  free-flight,  which  is  a  portion  of  an  elliptical 
ort,it  and  reentry,  which  is  the  portion  from  where  atmospheric  drag  becomes  sigraficant 
umil’missile  impact  The  powered-flight  phase  is  modeled  with  the  TBM  profile  polyno¬ 
mials  presemed  earlier.  The  ellipse  traced  during  free-flight  is  simulated  usmg  merUal 
two-body  mechanics  Atmospheric  drag  effects  during  the  r^emry  phase  are  not  spec.fi- 
cally  calculated  in  this  model,  but  are  somewhat  accounted  for  by  assuming  that  the  dis¬ 
tance  traveled  over  the  earth  from  re-emry  to  impaa  is  the  same  as  from  launch  to 
burnout  This  is  the  same  as  assuming  that  the  earth-central  angles  are  equal: 

e..=ei» 

The  speed  of  the  TBM  during  re-entry  is  assumed  to  be  unaffected,  however.  These 
apprtiximations  are  recognized  artificialities,  but  are  simpler  than  calculating  balhstic 
coefficients  of  various  TBM’s  and  modeling  atmospheric  density,  and  more  accurate  than 
assuming  the  missile  remains  on  its  elliptical  path  umil  impact.  Figure  15  illustrates  the 
geometry  and  some  of  the  quantities  used  in  the  equations. 
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Figure  15.  Ballistic  Trajectory  [Ref.  6] 

The  equations  contained  in  Chapter  6  of  Ref  6  concern  ballistic  trajectories  and 
are  based  upon  inertial  quantities  and  a  spherical  earth  model,  so  the  geocentric  position 
vector  and  inertial  speed  at  burnout  are  required.  Geocentric  latitude  is  generated  from 
(26): 


fbo  =  tan-’[(l-  f)^  tan(<t)fc„)] 


(26) 


The  earth’s  radius  at  burnout  is  calculated  with  (27): 

r  =  (27) 

V(1  -  0^  cos^  )  +  sin^  (<|>'fco ) 

These  values  are  then  transformed  into  earth-centered,  Carteaan  coordinates  of  the  TBM: 

cos(4>'^„ )  +  alt^  cos(<|>b^)]cos(A.^)  (28) 

Viu  =  cos(4>’b. )  +  alt^„  cos(<|)^)]sin(X^)  (29) 

Zt^  =  W  sin((l>^)  (30) 
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which  define  the  geocentric  TBM  position  at  burnout: 

*‘bo  ~  Ylk,  (^) 

In  an  inertial  reference  fi^e,  the  TBM  has  an  initial  eastward  velocity  at  launch 
due  to  the  earth’s  rotation,  Vo.  This  velocity  can  be  expressed  as: 

v„  =  rte«,^o>,cos(<|>;„)  (98) 

where  Oe  =  15®/hour,  the  rotation  rate  of  the  earth.  The  burnout  velocity  can  be  broken 
down  into  three  inertial  (UEN)  components.  In  the  inertial  frame,  the  initial  eastward 
velocity  is  added  to  the  eastward  component: 

V,  =V^sin(y^)  (99) 

cos(y^„)sin(a^)  +  v„  (100) 

Vk  =  Vb„cos(y,Jcos(a^)  (101) 

Inertial  quantities  are  denoted  with  a  bar,  “  ~  ”.  The  magnitude  of  this  vector  is  the  iner¬ 
tial  speed,  V^, : 

V,.  =  (102) 

The  inertia]  flight  path  angle  and  heading  may  now  be  calculated: 

Ybo  =  (103) 

_  ^  V  ^ 

ctbo  =tan-'j^^j  (104) 

Now  that  the  inertial  quantities  are  known,  a  non-dimensional  parameter,  Q,  is  defined  as 
the  squared  ratio  of  the  inertial  speed  of  the  TBM  to  circular  orbit  speed  at  that  position: 

V  *r 

Q  =  ■  -  ^  (105) 

M 

km^ 

where  p  =  398601.2  — the  gravitational  parameter  for  the  earth,  and  Tbo  =  |rbo|,  the 

sec 

magnitude  of  the  TBM  position  vector. 
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The  eccentricity  of  the  ballistic  “orbit”  is; 

e  =  Vl+Q(Q-2)cos'(ybo)  0^^) 

The  free-flight  “earth  central  angle”,  'T,  is  defined  by  Bate,  et  al.,  as  that  earth- 
central  angle  that  the  missile  traverses  between  burnout  and  re-entry  (see  Figure  1 5).  This 
definition  is  modified  by  adding  the  angle  traversed  during  reentry,  e„  to  the  BMW  de¬ 
fined  angle,  'F.  This  implements  the  assumption  that  the  TBM’s  tr^ectory  is  modified  by 
drag.  The  new  definition  for  the  angle,  T,  is: 


'F  =  2cos  ' 


l-Qcos^(yfc„)^ 


+  Q 


bo 


(107) 


The  eccentric  anomaly,  E,  and  the  semi-major  axis,  a,  of  the  ballistic  trajectory  are: 


E  =  cos 


e-cos^yj 

-  ^ 

1  -  ecoa  —  ,  , 
KlJ) 


(108) 


a  = 


'bo 


2-Q 


(109) 


Assuming  the  speed  of  the  TBM  is  unaffected  by  drag  during  re-entry,  the  time  of 
free  flight  (burnout  to  impact)  is  defined  as: 


tff  =  2  — [jt-E  +  esin(E)] 


(110) 


the  time  of  impact  fi-om  launch  is. 


tim  tbo  tff 


and  the  time  of  day  of  impact  is: 

Tim  ~  To  tim 

Using  the  Law  of  Corines  from  spherical  trigonometry,  latitude  at  impact,  4>'im,  is: 

4(;„  =  sin'’[rin(4);„)cos('F)  +  cos(27i-a^)] 

This  can  be  transformed  back  into  the  WGS-84  latitude: 


<|)i^  =  tan 


tan(4>L)^ 

(l-f)^J 


(111) 

(112) 

(113) 


(114) 
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This  method  is  not  an  exam  transformation,  but  it  is  a  dose  approximation  that  does  not 
affect  the  error  analysis,  A  better  model  should  be  used  if  more  precise  latitude  of  impact 

is  desired. 

Using  the  Law  of  Cosines  again,  and  after  some  algebra,  the  longitude  traversed, 


AX,  is; 


cos(y)-sin((|>l)sin(<|>U)] 

cos(4>;„,)cos(<|);,„)  J 


(115) 


So,  the  impact  lon^tude,  Xm,  is: 

=  (>>6) 

Better  models  for  the  ballistic  trajectory  could  possibly  be  employed.  The  goal  of 
this  analysis  is  not  to  precisely  determine  the  actual  impact  zone,  but  to  determine  the 
effects  of  error  sources  upon  the  result.  Thus,  the  effects  found  using  this  model  should 
be  applicable  to  other  ballistic  trajectory  models. 
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IV.  ERROR  SOURCES 


Noise  (error)  is  present  in  every  aspect  of  the  algorithm  presented  in  Chapter  3 
Every  measurement,  constant,  coefficient,  model,  equation,  approximation,  and 
assumption  impart  some  finite  quantity  of  uncettaimy  upon  the  result  Even  if  a  quantity  is 
exact  to  its  limit  of  accuracy,  the  simple  fact  that  it  is  represented  by  a  finite  numb^  of 
digits  limits  the  precision  of  that  quantity.  For  example,  if  a  ruler  has  one  sbctee 
gradations,  the  accuracy  of  any  measuremerns  made  with  the  ruler  is  ±one-thirty-second  of 
an  inch.  This  means  that  there  are  no  “perfect”  values  in  the  algorithm,  and  each  value 


used  is  a  source  of  uncertainty,  termed  an  error  source. 

Each  error  source  usually  has  one  or  more  underlying  causes  A  complete  error 
analysis  would  break  each  error  source  down  to  its  fundamental  level,  and  model  each 
level  correctly.  An  analogy  for  this  is  the  layers  in  an  onion;  The  onion  can  be  viewed 
externally  as  a  whole,  but  can  also  be  peeled,  layer  by  layer,  to  reveal  deeper,  underlying 
matter  that  support  the  external  skin.  In  this  section,  some  of  the  underlying  causes  of  the 
overall  error  sources  are  identified,  but  in  these  analyses,  only  the  overall  error  magnitudes 


will  be  modeled  and  analyzed. 

The  obvious  origin  of  any  errors  present  is  the  observational  data  since  it  is  the 
starting  point  for  the  algorithm.  These  data  are  physical  measurements  of  three  types: 
time,  focal  plane  measurements,  and  satellite  position  and  orientation. 

Time  errors  can  be  caused  simply  by  having  more  than  one  clock  being  referenced 

as  a  source  of  time  measurement  Imperfect  synchronization  between  clocks’  time  and 

time  passage  rate  are  obvious  error  sources.  Time  delays  caused  by  radio  transmission  of 
data  due  to  distance,  atmospheric  refraction,  and  relative  motion  Doppler  effects  may  add 
another  time  error  In  the  extreme,  the  source  of  time  error  can  be  traced  all  the  way  back 
to  our  measurement  of  time  passage  with  respect  to  the  celestial  sphere,  which  is  not 
fixed.  The  magnHude  of  time  errors  is  relatively  small,  and  getting  smaller  as  time 
measurement  improvements  are  being  made  continually. 
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The  largest  source  of  time  error  enters  the  algorithm  at  the  burnout  time  estimation 
phase.  The  10-second  sampling  rate  limits  the  accuracy  of  the  estimate  to  ±  5  seconds  for 
single  satellite  observations.  The  accuracy  of  the  profile’s  value  for  tnox  also  comes  into 
play  in  the  estimate.  The  burnout  time  is  crucial  to  determining  the  ballistic  trajectory, 
since  the  TBM  is  undergoing  its  greatest  acceleration  changes  during  the  final  seconds  of 
powered  flight.  Small  errors  in  the  estimate  of  burnout  time  thus  result  in  larger  errors  in 
the  estimate  of  impact  time  and  position. 

LOS  measurement  errors  can  arise  fi'om  many  noise  sources,  mainly  of  two  types: 
attitude  uncertainties  and  IR  radiation  measurement  errors.  The  attitude  of  the  spacecraft 
needs  to  be  precisely  known  (and  it  simplifies  the  process  if  it  is  very  stable).  A  28 
microradian  error  in  pitch  or  roll  pointing  accuracy  of  a  geosynchronous  satellite  equates 
to  a  1  kilometer  error  on  the  surface  of  the  earth,  at  least  35786  kilometers  away.  Any 
attitude  control  system  inaccuracy  effects  are  amplified  by  the  geosynchronous  altitude. 
Vibrational  noise  fi’om  spinning  momentum  wheels,  uncertdn  knowledge  of  the  mass 
properties  of  the  yaw-spinning  satellite,  star  catalog  and  star  sensor  measurement  errors, 
thruster  misalignments  and  disturbance  torques,  tachometer  errors,  and  mechanical 
misalignments  all  contribute  to  the  total  attitude  uncertainty.  The  knowledge  of  the 
telescope  boresight  axis  alignment  with  the  satellite’s  reference  fi’ame  is  defined  by  the 
design  and  manufacturing  of  the  DSP  satellite,  and  changes  slightly  with  thermal 
variations. 

The  IR  radiation  measurements  of  intensity  and  angle  of  arrival  (AOA)  also  have 
several  underlying  error  sources.  Locations  of  the  individual  photoelectric  cells  on  the 
focal  plane  array  are  recorded  in  what  is  termed  the  ‘Tocal  Plane  Vector  Table”  (FPVT). 
The  positions  of  the  PEC’s  change  as  the  satellite  heats  and  cools  with  varying  sun- 
satellite  orientations,  causing  a  warping  of  the  focal  plane,  but  the  FPVT  does  not  account 
for  the  changes  real-time.  In  addition,  detector  noise,  refiaction,  and  IR  attenuation  due 
to  clouds  and  water  vapor  all  add  to  the  total  measurement  error. 

Satellite  position  measurements  fi’om  AFSCN  are  used  to  update  orbit  ephemeris 
data.  The  measurements  are,  of  course,  inexact,  but  those  errors  are  compounded  over 
time  because  the  updates  occur  only  once  weekly.  The  ephemeris  data  is  propagated  to 
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determine  satellite  position  during  the  week  between  updates.  DSP  is  a  geosynchronous 
satellite,  so  hs  orbit  does  not  change  greatly  in  that  period,  but  errors  in  the  propagation 
model  can  cause  the  predicted  position  to  be  different  from  reality. 

Keeping  all  of  these  underlying  components  in  mind,  the  errors  are  modeled  in  the 
tactical  parameter  estimation  algorithm.  The  real-world  error  component  magnitudes,  bias 
and  random  distributions  may  be  different  from  those  simulated,  but  the  overall  effects 
manifest  themselves  in  a  marmer  close  to  the  error  models.  It  should  be  possible  to 
extrapolate  the  results  obtained  in  this  study  to  different  errors  by  properly  scaling  the 
different  magnitudes  and  distributions.  It  must  be  emphasized  that  since  the  goal  of  this 
study  is  not  to  exactly  determine  the  tactical  parameters  at  launch,  state  vector  at  burnout, 
or  impact  time  and  position.  The  purpose  is  to  determine  the  contribution  of  the  error 
sources  upon  these  values,  and  this  can  be  done  vdth  approximate  models  for  the  errors. 
If  the  intended  goal  is  to  calibrate  the  system  to  eliminate  bias  errors,  then  the  sources  of 
the  errors  must  be  determined  more  exactly  to  determine  what  is  observable  (measurable). 
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V.  NUMERICAL  ANALYSIS 


The  algorithm  described  in  Chapter  III  was  programmed  into  MATLAB™.  Three 
error  sources  (time,  satellite  position,  and  LOS)  were  simulated  and  added  to  the 
observational  data,  first  separately,  and  then  combined.  The  effects  of  the  errors  upon  the 
results  are  analyzed  at  three  points:  launch  position,  burnout  position,  and  impact 
position.  Five  Middle  Eastern  capital  cities  were  chosen  as  fictitious  launch  sites  to 
determine  the  effects  of  TBM  launch  position  upon  the  accuracy  of  the  results. 

The  model  for  time  error  is  a  random  uniform  distribution  between  0  and  1 
milliseconds.  For  satellite  position,  a  random  normal  distribution  with  zero  mean  and 
standard  deviation  of  200  meters  was  added  to  each  of  the  components  (R,  gha,  5).  This 
is  approximately  a  one  standard  deviation  sphere  with  radius  346  meters.  The  LOS  error 
was  modeled  as  a  random  normal  distribution  with  a  standard  deviation  of  5  microradians, 
added  to  both  P  and  ri  components.  Histograms  of  sample  error  distributions  are  shown 
in  the  following  three  figures. 
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Figure  16.  Histogram  of  Time  Error  -  Uniform  Distribution 
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Errors  added  to  the  observations  had  the  effect  of  creating  an  ellipse  of  incorrect 
data  points  around  the  true  portion,  temted  ».  “error  ellipse”  The  1^  ellipses  shown  in 
the  next  Bguies  are  determined  using  sutistical  methods.[Ref.  7]  For  the  launch  ellipse,  a 
matrix  of  two  columns  was  collected  during  the  simulation. 


lpts  = 


®«ror=0 


^0 _ OM.  ^0, 


'«ror999 


(117) 


where  Ipts  =  launch  points,  the  matrix  of  launch  longitudes,  Xo,  and  latitude,  (k  The  first 
row  is  the  generated  true  launch  position,  the  second  row  is  the  calculated  launch  position 
with  no  error  added  to  the  observational  data,  and  the  remaining  999  rows  are  the 
ealculated  launch  positions  with  errors  added  to  the  observational  data  The  6rst  two 
rows  are  eiror-free,  and  are  removed  to  produce  a  matrix  of  only  “corrupt”  positions  To 
determine  the  ellipse  dimensions,  a  and  b  (semi-major  and  semi-minor  axes  lengths),  the 
eigenvectors  and  values  of  the  covariance  of  this  matrix  were  calculated  Twice  the 
square  root  of  the  eigenvalue  diagonal  elements  produces  a  and  b. 


O  ic 

^eicenvalue(l,l)  ® 

IbJ 

=  2^ 

0  -^eigenvalue(2,2) 

and  the  eigenvector  is  the  direction  cosine  matrix  for  rotating  the  ellipses  from  the  primary 
axes.  The  impact  elUpses  are  calculated  similarly,  with  the  initial  data  matnx  being 
composed  of  impact  longitudes  and  latitudes,  instead  of  launch.  The  ellipses  can  then  be 


plotted  using  the  ellipse  equation; 


4+4=' 

a'  b' 


(119) 


and  letting  x  vaiy  between  ±a,  and  letting  y  be  the  dependent  variable.  Multiplying  the 
resuhant  (x,  y)  coordinates  by  the  eigaivector  matrix  rotates  them  to  the  coirect 
orientations  The  semi-axes  are  in  units  of  radians,  since  the  positional  values  are  in 
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radians.  Thus,  the  (x,  y)  plots  produced  are  the  longitude  and  latitude  of  the  points,  once 
the  values  are  converted  to  degrees. 

For  the  burnout  position,  the  same  method  is  used,  but  the  data  point  matrix  has 
three  columns,  the  third  being  the  altitude  of  the  IBM  at  burnout.  The  covariance  method 
works  as  well  in  three  dimensions,  producing  a,  b,  and  c,  the  three  semi-axes  for  the 
ellipsoid  formed.  The  coordinates  are  calculated  in  the  same  manner,  using  the  equation 
for  an  ellipsoid: 


2  2  2 

■;r+^+-j  =  i  (120) 

a  D  c 

and,  again,  rotating  by  pre-multiplying  by  the  eigenvector  matrix. 

The  following  figures  are  representative  of  all  the  simulated  cases,  but  this 
particular  case  is  a  300-km  TBM  launched  from  Baghdad  heading  0°,  with  the  combined 
erorrors. 


Launch  Position 


Figure  19.  Launch  Position  Error  Ellipse 

The  ellipse  drawn  on  top  of  the  data  points  encompasses  the  data  points  within  one 
standard  deviation  of  the  mean,  assuming  the  distribution  is  normal.  Since  the  time  error 
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has  a  uniform  distribution,  l-o  is  not  the  expected  66%  for  time  error  cases.  This  method 
is  used  anyway  for  lack  of  a  better  one  and  to  maintain  a  common  baseline  for 
comparison.  As  it  turns  out,  the  time  error  contribution  is  very  small,  and  thus  does  not 
effect  the  combined  error  case  distribution. 

The  position  at  burnout  is  plotted  in  three  dimensions,  using  altitude  as  the  third 
coordinate  to  form  an  ellipsoidal  error  volume. 

Position  at  BurnCut 


Figure  20.  Position  at  Burnout 

The  three  ellipses  show  the  outline  of  the  ellipsoidal  volume.  It  is  difficult  to 
perceive  the  ellipse  orientations,  but  they  are  mutually  orthogonal,  (necessarily,  since  they 
correspond  to  the  eigenvectors  of  distinct  eigenvalues). 
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The  impact  ellipse  is  similar  but  larger  than  the  launch  position  ellipse,  due  to  the 
propagation  of  velocity  errors  from  burnout  to  impact  and  the  error  in  the  burnout  time  . 

Impact  Position 


Figure  2 1 .  Impact  Position  Error  Ellipse 

The  three  previous  figures  illustrate  the  error  ellipses(oids),  but  showing  all  the 
data  obtained  in  this  format  is  unwieldy.  The  sizes  (areas  and  volumes)  of  the  positional 
errors  are  the  important  quantities  for  identifying  which  error  source  has  the  largest  effect. 
Thus,  the  data  have  been  processed  to  determine  these  quantities,  and  are  plotted  next. 

The  data  for  the  300-km  TBM  was  collated  by  city,  heading,  and  error  source. 
The  ellipse(oid)  dimensions  calculated  must  be  equated  to  the  actual  distance  on  the 
ground  to  determine  the  area  of  the  ellipses.  To  do  this,  the  angles  must  be  in  radians,  and 
are  multiplied  by  riocai-  The  distance  between  lines  of  longitude  shrinks  as  the  latitude 
moves  away  fi’om  the  equator,  so  the  longitudinal  distance  must  be  multiplied  by  the 
cosine  of  the  geocentric  latitude,  <!)'.  Thus  the  equation  for  ellipse  area  becomes: 

area  =  7tr,„^,,^ab  cos((t)')  (121) 
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The  ellipsoid  volume  is  not  as  simple,  so  a  parallelepiped  with  dimensions,  a  x  b  x 
c,  is  used  to  approximate  it: 

volume  =  abc  (122) 

These  areas  and  volumes  are  plotted  in  the  following  figures,  separated  by  error 
source  and  heading.  Each  launch  point  is  represented  by  a  different  color; 

Aden  =  blue 
Baghdad  =  black 
Damascus  =  green 
Riyadh  =  red 
Tehran  =  cyan 


Heading  (degrees) 

Figure  22.  Time  Error  Effects  Upon  Launch  Ellipse  Area 
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Figure  23.  Time  Error  Effects  Upon  Burnout  Ellipsoid  Volume 
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Figure  25.  Satellite  Position  Error  Effects  Upon  Launch  Ellipse  Area 


Figure  28.  LOS  Error  Effects  Upon  Launch  Ellipse  Area 
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Figure  30.  LOS  Error  Effects  Upon  Impact  Ellipse  Area 
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Burnout  Ellipsoid  Volume  (km"3)  Launch  Ellipse  Area  (km*2) 
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Figure  33.  Combined  Error  Effects  Upon  Impact  Error  Ellipse 
It  is  immediately  apparent  from  these  plots  that  the  launch  site  has  little  effect  upon 
the  size  of  the  error  ellipses  and  ellipsoids.  It  is  also  obvious  that  the  size  is  dependent 
upon  IBM  heading  and  observation  geometry.  The  relative  magnitudes  of  the  error 
effects  are  not  as  obvious,  so  a  plot  of  four  error  ellipses  is  shown  in  the  next  figure,  with 
each  ellipse  formed  by  a  different  error  source. 
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Launch  Position 


Longitude  (degrees) 

Figure  34.  Launch  Error  Ellipses  for  Each  Error  Source 
From  this  plot  it  is  easily  recognized  that  the  errors  can  be  ranked  in  order  of 

effect: 

1 .  focal  plane  errors 

2.  satellite  position  errors 

3.  time 

All  the  error  sources  together  produce  the  largest  ellipse,  as  would  be  expected  from  the 
superposition  principle. 


This  is  further  illustrated  by  Table  3.  Each  error  source  is  a  row:  “time”  is  the 
time  error  runs,  “satpos”  is  the  satellite  position  error  runs,  “LOS”  is  line  of  sight,  and 
“combin.”  is  the  combined  errors  runs. 


ERROR 

LAUNCH  ELLIPSE  AREA 
(km^) 

BURNOUT  ELLIPSOID 
VOLLnVIE  (km^) 

IMPACT  ELLIPSE  AREA 
(km^) 

min 

mean 

max 

min 

mean 

max 

min 

mean 

max 

time 

3.26e-8 

2.59e-7 

1.22e-6 

2.9e-II 

3.6e-10 

L33e-9 

1.57e-6 

2.6le-5 

BifHll 

3.77e-3 

5.27e-3 

8.68e-3 

6.04e-4 

I.I4e-3 

2.32e-3 

6.l2e-3 

7.92e-2 

LOS 

2.97e-2 

8.45e-2 

1.83e-l 

2.37e-2 

9.33e-2 

2.09e-l 

3.05e0 

1.02el 

1.95el 

combin. 

3.91e-2 

l.Ole-1 

2.97e-l 

4.44e-2 

I.25e-I 

4.50e-l 

3.35e0 

1.38el 

4.06el 

Table  3.  Error  Ellipse  Areas  and  Ellipsoid  Volumes 
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This  is  also  shown  graphically  in  the  following  plots; 


Minimum  Mean  Maximum 


Figure  35.  Launch  Ellipse  Area  Comparison 


Minimum  Mean  Maximum 


Figure  36.  Burnout  Ellipsoid  Volume  Comparison 


Minimum  Mean  Maximum 


Figure  37.  Impact  Ellipse  Area  Comparison 
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The  sum  of  the  error  sources’  positional  error  sizes  is  less  than  the  combined  error 
source’s  positional  error  size.  This  implies  a  nonlinear  interaction  of  error  effects  with 
each  other.  This  is  a  reasonable  result  considering  the  complexity  of  the  simulated  system, 
and  noting  that  most  real  systems  exhibit  nonlinear  behavior. 

Taking  advantage  of  the  ease  of  modifying  MATLAB™  code,  various  changes  to 
the  present  DSP  system  model  were  put  into  simulation  to  determine  the  effects  upon  the 
accuracy  of  the  results.  For  every  case,  Baghdad  is  the  launch  site  and  the  observational 
data  were  modified  by  combined  error  sources. 

The  first  case  is  the  control  case  with  nominal  system  parameters.  It  is  listed  as 
‘Bagcom”  to  represent  “Baghdad  combined  errors”,  and  is  used  as  a  baseline  case  for 
comparison.  The  second  case,  “Synchr.”  shows  the  effects  of  synchronizing  the  spins  of 
the  satellites  so  that  the  satellites  scan  a  single  area  of  interest  within  one  second  of  each 
other.  For  the  third  case,  “Molniya”,  the  satellite  at  70°  GHA  was  elevated  to  63.4°  decli¬ 
nation,  to  simulate  a  Molniya  +  geostationary  viewing  geometry.  The  remaining  cases 
were  simulating  faster  scan  rates  from  10  seconds  down  to  2.5  seconds.  The  numerical 


results  are  shown  in  Table  4. 


ERROR 

LAUNCH  ELLIPSE  AREA 
nan^) 

BURNOUT  ELLIPSOID 
V<XUME  (km^) 

IMPAC 

r  ELLIPSE 
(km^) 

.AREA 

min 

mean 

max 

min 

max 

max 

min 

mean 

max 

Bagcom 

4.21e-2 

l.Ole-l 

2.74e-l 

4.46e-2 

1.09e-l 

2.60e-l 

3.83e0 

1.40el 

4.04el 

Synchr. 

4.40e-2 

7.07e-2 

l.lOe-1 

2.09e-2 

4.89e-2 

8.61e-2 

3.94e0 

7.82e0 

1.52eO 

Molniya 

4.39e-2 

7.58e-2 

1.25e-l 

2.61e-2 

5.57e-2 

1.04e-l 

2.91e0 

7.96e0 

1.85el 

10  s  sc 

4.21e-2 

l.Ole-1 

2.74e-l 

4.46e-2 

1.09e-l 

2.606-1 

3.83e0 

1.40el 

4.04el 

7.5  s  sc 

2.62e-2 

7.44e-2 

1.95e-l 

2.53e-2 

4.55e-2 

1.09e-l 

4.03e0 

8.72e0 

2.08el 

5  s  sc 

2.48e-2 

5.02e-2 

1.52e-l 

L23e-2 

2.52e-2 

6.52e-2 

2.52e0 

5.84e0 

1.64el 

2.5  s  sc 

1.18e-2 

2.36e-2 

3.48e-2 

5.05e-3 

9.21e-3 

1.31e-2 

1.03e0 

2.95e0 

4.58e0 

Table  4.  Various  Case  Comparison 


The  effects  of  the  modifications  are  more  easily  observed  graphically; 
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Ellisoid  Volume  (km^) 


□  Baghdad  Combined 

□  Synchronized  Spins 

□  Moiniya 


Mean  Launch  Ellipse  Areas 


Figure  38.  Comparison  of  Mean  Launch  Ellipse  Areas 


Figure  39.  Comparison  of  Mean  Burnout  Ellipsoid  Volumes 
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Mean  Impact  Ellipse  Area 

Figure  40.  Comparison  of  Mean  Impact  Ellipse  Areas 

Both  modifications  have  the  effect  of  decreasing  the  areas  and  volume  by  about  a 
third.  This  is  the  expected  result  for  the  Molniya  case,  since  the  viewing  geometry  is  more 
three-dimensional  with  two  satellites  on  the  equator  and  the  middle  satellite  near  Molniyan 
apogee. 

The  synchronized  spin  results  were  a  surprise,  however.  The  effect  upon  the 
launch  ellipse  area  is  fine,  but  since  the  burnout  time  estimation  should  have  been  less 
accurate,  the  burnout  ellipsoid  volume  and  impact  ellipse  area  were  expected  to  remain 
unchanged  at  best.  Obtaining  simultaneous  observations  would  have  obvious  advantages, 
since  they  could  be  processed  differently  and  be  used  to  determine  the  position  exactly  for 
discrete  instants  in  time  during  the  boost  trajectory.  However,  these  observations  were 
not  processed  differently,  and  were  not  simultaneous  —  only  close  together  (within  one 
second  of  each  other).  Still,  the  result  may  be  a  manifestation  of  some  mathematical 
process  that  causes  the  increase  in  accuracy. 

The  following  plots  show  the  effects  of  increasing  the  scan  rate,  or  alternately, 
decreasing  the  time  between  observations.  This  allows  more  observations  to  be  made 
during  the  boost  phase. 
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JLllipsoid  volume  (km  ) 


0 10  s  scan 
M  7.5  s  scan 

□  5.0  s  scan 

□  2.5  s  scan 


Mean  Impact  Ellipse  Areas 

Figure  43.  Scan  Rate  Effects  on  Impact  Ellipse  Area 
The  effects  seem  to  be  linear  decreases  upon  the  launch  and  impact  ellipse  areas, 
and  quadratic  decreases  on  the  burnout  ellipsoids.  This  result  is  in  perfect  agreement  with 
intuitive  predictions. 

In  summary,  the  numerical  results  are  generally  in  accordance  with  expected  re¬ 
sults.  As  a  qualitative  check  of  the  quantitative  results,  the  next  chapter  is  devoted  to  a 
qualitative  analysis  of  a  simplified  case  for  comparison. 
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VI.  QUALITATIVE  ANALYSIS 


Starting  with  the  equations  presented  in  Chapter  3,  a  qualitative  analysis  has  been 
done  to  determine  the  expected  trends  in  the  numerical  analysis.  The  problem  is  simplified 
by  assuming  that  the  equations  are  exact  and  that  there  is  only  one  satellite.  The  eflFects  of 
time  errors  are  not  analyzed,  so  the  problem  is  to  determine  how  satellite  position  and 

focal  plane  error  effects  should  effect  the  results. 

TBM  geocentric  longitude  and  latitude  can  be  expressed  in  (XYZ)  coordinates; 


(}>'  =  tan"  ’ 


(12) 


These  (XYZ)  coordinates  can  be  expressed  in  terms  of  satellite  position  and  focal  plane 
coordinates  using  the  LOS  projection  equation  (8): 


r  = 


Yt 


=  R.+pe 


(8) 


Rewriting  the  terms  on  the  left-hand  side  in  terms  of  R,  gha,  5,  ti,  and  p  generates  a  rather 


lengthy  equation,  and  is  given  in  Appendix  B,  The  next  task  is  to  find  the  partial 
derivatives  of  the  TBM  position  with  respect  to  the  five  error  terms,  using  the  Cham  Rule. 

Taking  the  partial  with  respect  to  R,  for  example; 

8X  _  d\  dx  dX  dy  ^  dX  dz  (123) 

^~0xoR  6y5R  Sz^R 


d^’  dy  dz 

^”axaR^5yc?R  &0R 


(124) 


Again,  expression  this  equation  written  explicitly  in  terms  of  the  error  sources  becomes 
tedious,  and  can  be  found  in  Appendix  B.  The  results,  however,  are  interesting  and  are 
shown  as  plots.  The  horizontal  axis  pointing  left  is  r\,  the  horizontal  axis  pomtmg  right  is 
P,  and  the  vertical  axis  is  the  differential. 
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(125) 


The  partial  derivatives  with  respect  to  gha  turn  out  to  be  simple: 


(126) 


The  observations  one  can  make  from  analyzing  these  plots  is  that  most  of  the 
partial  derivative  plots  are  flat  and  close  to  zero.  When  riapproaches  8.6°,  however,  some 
of  the  values  begin  to  get  large.  This  happens  when  the  observed  IR  event  is  close  to  the 
limb  of  the  earth  as  viewed  from  a  geosynchronous  satellite.  There  also  appear  to  be 
transitions  when  t|  approaches  8.6°  and  b  is  0°  or  180°. 

Figures  44  and  45  show  that  satellite  radius  error  effects  are  of  such  small 
magnitude  that  they  are  neglible.  Figures  46  and  47  show  that  satellite  declination  error 
effects  are  more  important  to  determine  latitude  than  longitude,  except  when  r| 
approaches  8.6°.  Figures  48  and  49  show  that  P  LOS  errors  affect  latitude  determination 
more  than  longitude,  except,  again,  when  r|  approaches  8.6°.  Figures  50  and  51  show  that 
T|  LOS  errors  have  the  largest  effect  of  all,  and  the  partials  exhibit  the  same  behavior  near 
the  limb. 

These  plots  of  formulae  give  an  indication  of  what  magnitude  the  error  ellipse 
should  have,  given  the  magnitude  of  the  error  is  known,  by  using  that  magnitude  in 
equations  similar  to  (127): 


,,  ax  fax  ax  ax  ay  ax  azV„ 

aR  V  ax  aR  ^  aR  az  aRy 


(127) 


where  R  can  be  any  of  the  five  error  quantites  (R,  6,  gha,  t|,  or  P). 

To  test  this  assumption,  substitute  the  magnitude  of  the  errors  used  in  the 
simulation,  and  compare  the  result  with  the  numerical  analysis.  The  errors  are: 

AR  =  200m  (128) 


A5  =  tan"' 


0.200  km 
42164.17  km 


(129) 
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Agha  =  tan"' 


0.200  km 
42164.17  km 


(130) 


Ati  =  5  ^radians  (131) 

Ap  =  5pradians  (132) 

The  generalized  equations  for  the  partial  derivatives  of  the  geocentric  longitude 
and  latitude  can  be  used  to  find  the  derivatives  for  a  specific  location  by  using  specific 
focal  plane  coordinates.  For  Baghdad,  viewed  from  the  satellite  at  70°  gha,  in  particular, 
these  coordinates  would  be; 

T]  =  0.1216  radians  (133) 

B  =  3.8546  radians  (134) 


Ti  =  0.1216  radians 
P  =  3.8546  radians 


Using  these  values,  the  partial  derivatives  with  respect  to  the  five  error  sources 


become: 


a 

rad 

-0.573  — 

(135) 

ap~' 

rad 

d^'  _ 

rad 

-0.753 

(136) 

ap 

rad 

dX  _ 

rad 

-13.788  — 

(137) 

dr] 

rad 

at)'  _ 

.r..  rad 

5.96 - 

(138) 

an 

rad 

dX  _ 

rad 

0.656  — 

(139) 

a5 ' 

rad 

aij)'  _ 

„  rad 

-0.658  — 

(140) 

as 

rad 

_a  _ 

-3.663  xl  0" 

*  rad 

(141) 

aR 

m 

:  1.583  xl  0"* 

rad 

(142) 

aR 

m 

dX 

11 

II 

(143) 

agha 

rad 
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(144) 


-^  =  0  — 
dgha  rad 

Multiplying  (135)  through  (144)  by  the  appropriate  error  quantities  produces  the 


geocentric  longitude  and  latitude  error  values. 

=3.166x10“*  rad  (145) 

=  0  rad  (146) 

A<|>;  =-3.121x10"*  rad  (147) 

A<|);  =2.98x10"*  rad  (148) 

A<t)' =-3.765x10“*  rad  (149) 

AX  R  =-7.326x10“*  rad  (150) 

AXgt.  =4.743x10^  rad  (151) 

AXg  =-3.112x10^  rad  (152) 

AX ^  =-6.894x10"*  rad  (153) 

AX p  =-2.865x10“*  rad  (154) 


Grouping  these  in  terms  of  satellite  position  and  focal  plane  and  summing  the  absolute 


values: 

+4d>i,.  +A*;  =6.287x10^  rad  (155) 

=  AX,  +  AX,,.  +  AX,  =  1.518  X 10-’  rad  (156) 

Ad-L.,,,..,™  =  Ad;  +  Ad.;  =  3.357  x  lO''  rad  (157) 

AX,..,,,..™  =  AX, +4X,  =  7.181x10-’  rad  (158) 

These  values  are  in  radians,  Avhich  can  be  equated  to  kilometers  on  the  surface  of 
the  earth  by  muhiplying  latitude  by  the  earth’s  radius  and  longitude  by  the  earth’s  radius 
and  the  cosine  of  the  latitude.  For  this  analysis,  earth’s  radius  is  assumed  to  be  6378.137 
kilometers. 

Ad.;..,.,.,,.™  =(6.287x10'’'  rad) 6378.137  km  =  4.010x10-’ km  (159) 

AX™,.,™.™  =(1.518x10-’  rad)cos(33.333>378.137  km  =  8.089x10-’ km  (160) 
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A4»U.pUa.«-  =  (3-357  X 10-’  rad)  6378.137  km  =  2.141  x  lO'*  km  (161) 

=  (7-181  X  lO'*  rad)  cos(33.333“)6378.137km  =  3.827  x  10''  km  (162) 

Assuming  that  these  values  are  equivalent  to  the  semi-axes  of  an  error  ellipse,  multiplying 
longitude  error  by  latitude  error  and  by  tc  gives  the  area  of  the  error  ellipses; 

ellipse  area^ui.p„i^«^  =  1.019  x  10'^  km^  (163) 

ellipse  areaf^,p^.„„  =  2.574  x  10  '  km^  (164) 

The  focal  plane  error  effects  are  25  times  greater  than  the  eflFects  of  the  satellite  position 
error.  The  relative  magnitudes  agree  with  the  numerical  results,  but  the  specific  values 
obtained  here  are  greater  than  those  obtjuned  numerically  by  a  factor  of  about  three.  The 
numerically  obtained  values  were: 

launch  ellipse  mean  area„„Bitepo.i^^„  =  5.27  x  10“*  km^  (165) 

launch  ellipse  mean  area^^, puoeenor  =  8.45  x  10'^  km^  (166) 

This  difference  is  expected,  since  this  analysis  does  not  account  for  stereo  viewing,  the 
least  squares  iteration  process,  or  time  effects. 

This  result  shows  that  it  is  possible  to  compute  error  ellipse  areas  using  a 
simplified  qualitative  analysis.  The  relative  order  of  the  magmtude  of  the  effects  is  in 
agreement  with  predictions,  and  verifies  the  results  obtained  numerically.  To  do  the  same 
analysis  for  burnout  volumes  and  impact  areas  would  require  the  addition  of  time  errors, 
and  trajectory  equations.  This  increases  the  difficulty  of  this  “back  of  the  envelope” 
analysis,  and  is  not  treated  here. 
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VII.  CONCLUSION 


The  TALON  SHIELD  /  ALERT  state  vector  estimation  algorithm  was  correctly 
modeled  in  MATLAB’^.  System  errors  were  modeled  and  introduced  into  the  algorithm 
to  determine  the  eflFects  upon  the  accuracy  of  the  final  results.  A  simplified  qualitative 

analysis  verified  the  quantitative  results. 

The  three  questions  posed  in  the  first  chapter  are  now  answered: 

1 .  The  errors  in  the  system  are  broadly  categorized  as  errors  in  time,  satellite 
position,  and  LOS. 

2.  The  relative  magnitudes  of  the  error  effects  listed  largest  to  least. 

A.  LOS  -  using  a  zero-mean,  5  pradian  standard  deviation  normal  error 
distribution  in  focal  plane  coordinates 

B.  satellite  position  -  using  a  zero-mean,  200  meter  standard  deviation  normal 
error  distribution  in  satellite  positon  coordinates 

C.  time  -  using  a  uniform  error  distribution  between  zero  and  1  millisecond. 

3  The  effects  seem  to  behave  independently,  since  the  principle  of  superposition 
works.  The  effects  of  the  combined  errors  are  slightly  greater  than  the  sum  of 
the  individual  effects. 

Additional  runs  were  made  to  determine  the  effects  of  various  geometries,  spin  rates,  and 
observation  synchronization. 

There  are  many  tasks  left  to  accomplish  in  this  area  of  research.  The  MATLAB 
code  can  be  optimized,  additional  parameters  and  system  changes  can  be  simulated,  and 
more  detailed  analysis  of  the  results  can  be  made.  The  error  sources  can  be  modeled  more 
exactly  by  analytically  “peeling  the  error  onion”,  or  by  using  actual  DSP  values  and 
measurements.  The  results  obtained  could  be  further  analyzed,  possibly  to 

Simulating  the  DSP  system  with  computer  code  has  proven  to  be  an  effective  tool 
for  error  analysis,  which  is  a  necessary  first  step  for  determimng  how  best  to  improve  the 
present  TBM  detection  system  or  modify  the  design  of  a  follow-on  system. 
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APPENDIX  A.  “A”  MATRIX  DERIVATION 


The  A  matrix  is  made  up  of  partial  derivatives  of  focal  plane  coordinates  with 
respect  to  the  tactical  parameters.  These  elements  can  be  determined  qualitatively.  Exact 
equations  are  not  required  for  the  iteration  process  to  converge,  so  small-angle 
assumptions  are  made  to  simplify  the  resulting  equations.  The  benefit  of  using  these 
approximations  for  the  matrix  elements  is  the  ease  of  computation  and  thus  a  reduction  of 

computing  time. 

The  A  matrix  is  divided  into  three  components,  Ai,  A2,  and  A3.  The  elements  of 
Ai,  and  A:  will  be  derived  in  order,  using  MATHCAD  5.0  Plus.  The  A3  matrix  is  derived 
in  Chapter  3.  The  elements  of  the  first  column  of  Ai  are  partials  of  latitude,  <|). 

^  =  ^  -  acosicos(0)  sin\4>  oj  sin(6)  coS;^^  Q^-cos^a  Q  j  j 


Taking  the  partial  derivatives; 

■  sin(0)  sin  (b  q  ^  cos(0)  cos:^>  Qi  cos,a  q; 


d0  i  .  2 

_|1  ^cos(0)  sin;(h  q1  ^  sin(0)  cos^^  O/  ' 

Using  the  small  angle  assumption  that  sin(9)=0  and  cos(0)=l 

d*  cos.  ♦  0; 

d9  , - 2  ' 

.11 

^■efi 


dd  ref! 

d=(l  -  1.5  L)  dp 

dp*ao  + Tq]  +  a2(^T- Toj  +a3iT-TQj  +a4(^T-TQj 
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dTo 

dd 


a(l-  1.5L)[-ai -2a2  j^T- Tq)- 3  a3  i^T-To) 

1.5-[ao  +  ai-(T-Toi  +  a2-(T-To/  +  a3-^T-Tof +  a4-iT-To)  J 


dL 


d(t)  d^  d9  dd  dd 

H^"d9dddTo  dTo 

d^_d4  ^  dd  -'l  ^  dp  cos^g  q: 
dL  d9  dd  dL  r^g- 


d0  0 


d0 


aO 


d>-0 


dho 


d<t»  n  • 

- a  6  sin  ,  a  0 

da  0 


The  second  column  is  composed  of  partials  of  longitude,  X,: 


,  sin(0)  sin.a  q 
Xa).  Q+asin'  — 


COS(0  ) 


d9 


_ 

2  sin;  a  of 
1  sin(e)"— —  2 

COS(0) 


cos(6) 


sin  a  o; 
cos(0 ) 


sin  ,  a  0 , 
cos(0 ) 


A  ^  dd  0.)  dd  _  "'"v^  0)  dd 

dT^'dedddTo  r  jjf cos(» )  dTo  r  J(rcos^♦  qj  dTg 


Assuming  ^  is  approximately  <l>o: 
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dl  dl  dd  dd  -1  Sdpsmi^aoj  -l.Sdpsin^ 


dL  de  dd  dL  rg^cos((^)  rggpcos^*  q; 


dl  _ 
d^  0 


•sin(e) 


^1“  Oj  .  .  sin(0) 


1  -  sin(6) 


2  o; 

cos(^)^ 


COS((ll  ) 


sin(4)  -^ssiii(e)  sin^a  q  , 
d^  0 


COS((|)  ) 


dl 


dlo 
dl 


-1 


dh 


=0 


0 


d/w 

da  0 


I  2  0  / 

1  sin(e)  : 


•sin{9) 


cosiart;  BcoS'ao) 

\  ^ \ 

COS((t))  cos  oj 


cos(^ ) 


And  the  third  column  elements  are  the  partials  of  height,  h; 

h=h  Q  ‘  ( 1  -I  L)  h  p 


b3T  To)+b4iT  Tq, 


1  D- 1  -  b  ]  -  2  b  2- ;  T  -  T  0 )  -  3  b  3-  i  T  -  T  ^  -  4  b  4-  ( T  -  T  0 


dT 


dh 


=hi 


dL 

dh 

d(i)  0 

dh 

dh 


=0 


dh 


0 
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The  A2  matrix  is  the  partials  of  (XYZ)  with  respect  to  (<J>,  X,  h); 

Xt,  =  [rio«u  cos(<t»\  )  +  alt^  cos(4>j]cos(>.^) 

yT,  =  [•’loe.ik  cos(<t)'J  +  alt^  cos((|)  j]sin{X^) 

Zt,  =  sin(4»\  )  +  alt^  sin((|>  J 

Approximating  pocai  with  reff ,  alt  with  h,  (f>'  with  <[),  and  talcing  the  partials: 

■^  = -(r^  ^  h)sin((t»)cos(?i) 
oj) 

=  -(feff  +h)c05((|>)sin(?i) 

CA 

^  =  cos(<t))cos(X) 

an 

=  -{r^ff  *  h)sin(<l))sin(A) 

^  +  h)cos((t>)cos(?i) 

CA 

dy 

—  =  cos(<j))sin(}i) 


^  =  (r,^+h)cos(<|») 

^=0 

dk 


=  sin(<|)) 
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APPENDIX  B.  QUALITATIVE  ANALYSIS  EQUATIONS  AND  PLOTS 
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Z2 
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Differentiating  with  respect  to  R: 


t^8 


The  geocentric  longitude  and  latitude  can  be  expressed  in  terms  of  (XYZ) 
coordinates: 


iy  ^ 

JlMatan  i  - 1 

U/ 


dx 

’"'Hi. 

dX_ 

1 

dy 

L  \  ^  /  J 

4  'saUn  j 


W^  +  y^ 


+ y^l 


;  2  2,  i 

\X  +y  ■  j 


?!  2  ! 
I-  ,  i  ,  7.  I 


dz”  ! 


ix'-y".' 


The  differentials  of  the  geocentric  longitude  and  latitude  can  be  determined  using 

the  Chain  Rule.  Using  longitude  and  satellite  radius  as  example  quantities: 

d>,  _  dX  dx  dX  dy 
dR  dx  dR  dy  dR 

Substituting  the  following  quantities  into  the  previous  equations: 

R  =42164.17kii] 

=70  deg 

8  =0 

‘‘local  '^^^* 
r\  =  varies  0°  to  8.5° 
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p  =  varies  0°  to  360° 


This  gives  expressions  for  the  variations  of  the  geocentric  longitude  and  latitude 
due  to  errors  in  LOS  projections  from  a  DSP  satellite  with  a  Greenwich  hour  angle  of  70°. 
These  expressions  are  transformed  into  plots  for  examination. 

The  following  plots  show  the  differentials  of  geocentric  longitude  and  latitude  with 
respect  to  the  five  error  sources.  The  horizontal  axis  pointing  to  the  left  is  i],  the 
horizontal  axis  pointing  toward  the  right  is  3,  and  the  vertical  axis  is  the  differential 
evaluated  at  each  (t|,  ft)  coordinate. 
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-t).998 


•dy.  ^ydP:  : 


31.746 


-31.746 


241.145 


“241.145 


I  fill 
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APPENDIX  C.  MATLABCODE 


The  following  is  the  MATLAB’’^  program  used  to  simulate  the  TALON  SHIELD 

/  ALERT  state  vector  estimation  algorithm.  It  is  commented  in  blue  font. 

%  .  Martin  Beaulieu 

%  DSP  Estimation  Algorithm 

%  Initialize  variables 


clear; 

format  long; 

numpoints=1000; 

reff^6371; 

f=l/298.257; 

re=6378.137; 

a0=0.0749291380897402; 

al  =-0.0501 64742464279 1 ; 

a2=0.0041947951949387; 

a3=-0.0000141392516462238; 

a4=8.51293969732759E-07; 

b0=0.884282320240989, 

bl=-0.129163814041924; 

b2=0.0137288843548327; 

b3=-0.0001 59307288768991 , 

b4=  1.3693  8557528 199E-06; 

tmax=62.5; 

mu=398601.2; 

d2r=pi/180; 

r2d=l  80/pi; 

we=15*d2r/3600, 


%  clear  memor>’ 

%  use  1 5-digit  numerical  format 
%  number  of  data  points  to  run 
%  spherical  earth  radius  (hm) 

%  \VGS-84  flattening  factor  for  oblate  earth 
%  oblate  earth  equatorial  radius  (km) 

1 0  coeflicients  in  range  profile  equation 


“'0  coeflicients  in  height  profile  equation 


TBM  maximum  motor  burn  time  (seconds) 
“o  earth  gravitational  constant  (km '  3/sec '2) 
%  degrees  to  radians 
%  radians  to  degrees 
%  earth  rotation  rate  (15  degrees'hour) 


%- 


%  This  section  generates  the  obsersation  table  (time,  az.  el,  gha,  dec,  and  radius).  The  first 
%  step  is  to  choose  the  tactical  parameters,  then  calculate  how  that  missile  profile  would 
%  appear  to  an  orbiting  sensor. 

%  Initial  tactical  parameters: 

for  ii=0;  1 1  %  loop  through  1 2  headings  around  compass  rose 


T0=100; 

%  launch  time  of  day  in  seconds 

L=0; 

%  loft  parameter 

%latO=  1 2 . 783  *d2r;lon0=45 .050*d2r; 

“/o  .Aden,  Yemen 
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Iat0=33.333*d2r,lon0=44.433*d2r; 
%latO=33 .500*d2r,lon0=36.3 19*d2r; 
%lat0=24.633*d2r;lon0=46.717*d2r, 
%latO=3  5 .666*d2r;lon0=5 1 .433  *d2r; 
h0=0; 

hdg0=ii*pi/6; 

data=[TO  L  lonO  latO  hO  hdgO  zeros(l,9)]; 

%  Generate  observation  time  matrix: 

obsl  =T0+30+5  *randn(  1 ); 
dt  1 = lO^randC  1  );dt2=  1 0*rand(  1 ); 
ttime=[obsl  ;obs  1  +dt  1  ;obs  1  +dt2] ; 
ttime=sort(ttime); 
nextobs=obsl+10; 

j=i; 

while  nextobs  <=  tmax+TO, 
ttime=[ttime;ttimeO)+ 1 0]; 

j=j+U 

nextobs=ttime(j)+l  0; 
end 

n=length(ttime); 

“o  Generate  saiellite  position  matrix 

gha=[10*d2r;70*d2r,105*d2r]; 
satord=ceil(6*rand(  1 )); 
if  satord  =  1 

tgha=[gha(l  ),gha(2);gha(3)], 
elseif  satord  —  2 
tgha=[gha(l);gha(3);gha(2)3; 
elseif  satord  =  3 
tgha=[gha(2);gha(  1  );gha(3)] ; 
elseif  satord  ==  4 
tgha=[gha(2);gha(3);gha(  1 )] , 
elseif  satord  =  5 
tgha=[gha(3),gha(  1  ),gha(2)] ; 
elseif  satord  ==  6 
tgha=[gha(3);gha(2);gha(l)]; 
end 

tdec=[0;0;0]; 

tradius=[421 64. 1 7;42 1 64. 1 7;42 1 64. 1 7]; 
for  i=4;n 

tgha=[tgha;tgha(i-3)]; 
tdec=[tdec;tdec(i-3 )] ; 


%  Baghdad-  Iraq 
%  Damascus,  Syria 
%  Riyadh.  Saudi  .Arabia 
%  Tehran.  Iran 
%  launch  altitude 

%  launch  heading  in  30  degree  steps 
%  true  tactical  parameters 


%  first  observation  time  -  clouds,  vapor,  etc 
%  two  random  time  intervals,  0<dt<]0  sec 
%  first  3  observation  times 
%  put  times  in  chronological  order 
%  next  sequential  obserx'ation 
%  counting  index 

%  if  next  obserx  ation  is  during  motor  burn 
®  o  add  next  obsen  ation  to  time  matrix 
increment  counter 
®  o  next  sequential  obser\  ation 
®  0  loop  until  matrix  is  completed 
°o  number  of  observ  ations 


“fi  (ireenwich  hour  angle 
0  randomh  select  the  order  in 
®  0  which  satellites  obsen  e  TBM 


%  declination 

%  geosynchronous  radius 

%  make  satpos  matrix  the  same  size  as  time 
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tradius=[tradius;tradius(i-3,:)]; 

end 

satpos=[tradius.  *cos(tdec).  *cos(tgha) 
tradius.  ’•'cos(tdec),  *sin(tgha) 

tradius.*sm(tdec)];  %  satellite  position 


%  Apply  profile  to  observation  times: 


t=ttime-T0; 

dp=aO+al  *t+a2*t.'^2+a3  *t/"3+a4*t.M; 
d=(l-1.5*L)*dp; 

hp=bO+bl  •t+b2*t/'2+b3  *t.^3+b4*t.M; 
h=(l+Lrhp; 


%  time  of  flight 
%  nominal  profile  distance 
%  loft-modified  distance 
%  nominal  profile  height 
%  loft-modified  height 


%  Generate  target  position: 

theta=d/reff,  earth-central  angle 

lati=pi/2-acos(cos(theta)  *  sin(latO)+sin(theta)*cos(latO)*  cos(hdgO));®  o  latitude 

long=lonO+asin(sin(theta)*sin(hdgO)./cos(lati)),  %  longitude 

alt=hO+h;  altitude 

gcl=atan((l-f)^2*tan(lati));  geocentric  latitude 

rloc=re*(l-f)  /sqrt((l-f)^2*cos(gcl).^2+sin(gcl)."^2);  ®/o  local  radius 

tgtpos=[(rloc.  *cos(gcl)+alt.  *cos(lati)).  *cos(long) 

(rloc.  *cos(gcl)+alt .  *cos(lati)).  *  sin(long) 

rloc.*sin(gcl)+alt.*sin(lati)];  %  target  position 


°o  Generate  line-of-sight  vector  and  transform  into  focal  plane  coordinates 

delta=tgtpos-satpos;  “o  line-of-sight  vector 

r3=[];r3t=[];UEN=[];  %  reset  variables 

for  i=l:n 

rotl=[cos(tgha(i))  -sin(tgha(i))  0;sin(tgha(i))  cos(tgha(i))  0;0  0  1]; 

rot2=[cos{tdec(i))  0  -sin(tdec(i));0  1  0;sin(tdec(i))  0  cos(tdec(i))]; 

rot3=rotl  *rot2; 

r3t=[r3t;rot3']; 

r3=[r3;rot3]; 

UEN(i,:)=(r3t(3*i-2:3*i,:)*delta(i,:)')';  ®/b  Up-East-North  coordinates 

end 

tbeta=rem(-pi/2-atan2(UEN(:,3),UEN(:,2))+(2*pi),(2*pi));%  true  azimuth 
teta=atan(sqrt(UEN(:,2).'^2+UEN(:,3).''2)./(-UEN(:,l)));  %  true  elevation 


°'o  Using  generated  data  (ttime,  theta,  teta.  tgha,  tdec,  tradius), 

® b  compute  tactical  parameters  using  least-squares  iteration  method: 
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for  index=l  ;numpoints; 
index=round(index); 


%  ensure  index  is  an  integer 


%  Initialize  variables: 

clear  time  radius  gha  dec  beta  eta  TO  L  latO  lonO  hO  altO  hdgO  t^ 
sse=l ;dsse=l ;ssel=l ;failindex=0;r3t=[];r3=[];  %  reset  vanables 

%“adiuserror 

.isa.e,,,.eg.,  error 

j  -«r«crn  iv  %  satellite  dec  eiTor 

errdec=zeros(n,l),  *1, 

,  ^  /«  1  %  azimuth  error 

erTbeta=zeros(n,l),  "  7 

.  1 1-  /o  elevation  error 

erreta=zeros(n,  1 ), 

°rb  Chanue  time,  beta,  eta,  and  satpos  by  some  error,  epsilon 


if  index  >  1 

errtime^  1  e"3  *(rand(n,  1  )"0. 5), 

errrad=0. 2*randn(3 , 1 ) ; 
errgha=atan(0. 2/421 64. 1 7)*randn(3 , 1 ), 
errdec=atan(0.2/42164. 1 7)*randn(3,l), 
for  i=4;n 

errrad=[errrad;errrad(i-3)]; 
errgha=[errgha;errgha(i-3)]; 
errdec=[errdec;errdec(i-3 )] ; 

end 

errbeta=5e-6*randn(n,  1 ); 
erreta=5e-6*randn(n,  1 ); 
end 

time=ttime+errtime; 

radius=tradius+errrad, 

gha=tgha+errgha; 

dec=tdec+errdec; 

beta=tbeta+errbeta; 

eta=teta+erreta; 


o-o  <=-.0.001  second 
%  variance  of +-200  meters 
®  o  variance  of  “^-200  meters 
“o  variance  of  ~-200  meters 
» 0  fix  satellite  positions 
during  single  run 


“  o  \  ariance=  --5microradians 
®  o  \'ariance=  -r--5microradians 

“  0  time  =  truth  -  error 
°  0  radius  =  truth  "T  error 
gha  =  truth  -  error 
rio  dec  =  truth  ■+  error 
°/o  beta  =  truth  ^  error 
eta  =  tnjth  -  en  or 


%  Compute  target  position: 

satpos=[radius.  *cos(dec).  *cos(gha) 
radius.*cos(dec).*sin(gha) 

radius.*sta(dec)];  “■'»  =^atellrte  position 

for  i=l  n  rotation  matrices 

rotl=[Ms(gha(i))  -sm(gha(i))  0;sm(gha(i))  cos(gha(i))  0;0  0  1], 
rot2=[cos(dec(i))  0  -sin(dec(i));0  1  0;sin(dec(i))  0  cos(deo(i))]; 
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rot3=rotl*rot2; 
r3t=[r3t;rot3']; 
r3=[r3;rot3]; 
uen=[-C50s(eta(i)); 

-sin(beta(i))*sin(eta(i)); 
-cos(beta(i))*sin(eta(i))] , 
ehat=T0t3*uen; 

chi=pi-asm((norm(satpos(i,:))*sin(eta)/reff)); 

rho=TeflE7sin(eta)*sin(eta+chi); 
los=Tho*ehat; 
tgtpos=los+satpos(i,;)'; 
lat(i)=atan(tgtpos(3)/sqrt(tgtpos(l)^2+tgtpos(2)^2)); 
lon(i>=atan2(tgtpos(2),tgtpos(  1 )); 
end 


%  Up-East  -North  coordinates 
%  line  of  sight  direction 
%  computation  of  length 
%  for  los  vector 
%  iine-of-sight  vector 
%  target  position 
%  latitude 
®'o  loncitude 


°o  Intial  estimate  of  tactical  parameters 

latO=(lat(  1  )+lat(2)+lat(3))/3; 
lonO=(lon(  1  )+lon(2)+lon(3))/3 ; 
latbo=(lat(n-2)+lat(n- 1  )+lat(n))/3 , 
lonbo=(lon(n-2)+lon(n- 1  )+lon(n))/3 , 
hdgO=rem(pi/2-atan2(latbo-latO,(lonbo-lonO) 
T0=time(l)-20, 

L=0; 

h0=0; 

tp=[TO;L;latOdonO;hOdtdgO], 


°  o  first  obs  latitude 
%  first  obs  longitude 
°'o  latitude  at  last  obs 
%  longitude  at  last  obs 

*cos(lat0))+(2*pi),(2*pi)),  ®  olnch  lieading 
%  launch  time 
“'0  loft  parameter 

launch  height  above  WGS-84  ellipsoid 
°  0  tactical  parameters  matrix 


%  Iterate  on  initial  estimates  until  the  difference  between  the  sum  of  the 
%  squares  of  the  errors  between  two  consecutive  runs  is  <=10*eps: 


while  dsse  >  10*eps 
duv=[];A=[]; 

T0=tp(  1  );L=tp(2);latO=tp(3); 

Ion0=tp(4);h0=tp(5);hdg0=rem(tp(6),2*pi); 

t=time-T0; 

dp=aO+al  *t+a2*t.^2+a3*t.^3+a4*t.''4; 
d=(l-1.5*L)*dp; 

hp=bO+b  1  *t+b2*t.^2+b3  *t.''3+b4*t.M, 

h=(l+L)*hp; 

theta=d/reff; 


%  loop  until  dsse  <  10*eps 
“'0  reset  variables 
%  update  tp  matrix  values 
with  dtp's  added 
%  time-of-flight 
%  nominal  profile  distance 
%  lofi-modified  distance 
%  nominal  profile  height 
%  loft-modified  heiaht 
%  earth-central  angle 


lati=pi/2-acos(cos(theta)*sin(latO)+sin(theta)*cos(latO)*cos(hdgO));'?o  geodtic  latitude 
long=lonO+asin(sin(theta)*sin(hdgO)./cos(Iati));  %  longitude 
alt=hO+h;  %  altitude 

gcl=atan((  1  -f)^2*tan(lati));  °o  geocentric  latitude 

rloc=re*(l-f)./sqrt((l-f)^2*cos(gcl)/'2+sin(gcl)/"2);°o  local  radius 
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tgtpos=[(rloc.  *C5os(gcl)+alt.  *cos(lati)).  *cos(long) 
(rloc.  *cos(gcl)+alt.  *cos(lati)).  *sin(long) 


rloc.  *sin(gcl)+alt.  *sin(lati)]; 
delta=tgtpos-satpos; 
fori=l:n 

UEN(i,:Kr3t(3*i-2:3*i,:)*delta(i,:)')'; 

end 

uhat=-UEN(;,2)./UEN(:,l); 

vhat=-UEN(:,3)./UEN(;,l); 

u=-tan(eta).  *sin(beta); 

v=-tan(eta).  *cos(beta); 

du=u-uhat; 

dv=v-vhat; 


%  target  position 
%  line-of-sight  vector 
%  rotate  into  UEN  coordinates 


%  estimated  focal 
%  plane  coordinates 
%  actual  focal 

%  plane  coordinates 

%  difference  between 

%  estimate  and  actual 


’o  Compute  A  matrix: 


cl=cos(hdgO)/reff,c2=sin(hdgO)/reflE/cos(latO);  constants  in  A  matrix 
c3=sin(hdgO)/reff,c4=cos(hdgO)/reflD'cos(latO), 

dddT0=-(l-1.5’''L)*(al+2*a2*t+3*a3*t.^2+4*a4*t.'^3);  -horizontal  speed 
dhdT0=-(l+L)*(bl+2*b2*t+3*b3*t.^2+4*b4*t.''3);  ®o  -vertical  speed 
dudUEN=[UEN(:,2)./UEN(:,l)/'2  -l./UEN(:,l)  zeros(size(time))]; 
dvdUEN=I^N(;,3)./UEN(;,l)/'2  zeros(si2e(time))  -l./UEN(:,l)]; 

%  change  in  focal  plane  coordinates  ^^rt  changes  in  L  EN  coordinates 
for  i=l:n  ^  o  construct  A  matrix 

duv=[duv,du(i),dv(i)];  "  o  error  v  alues  matrix 

Al=  [dddT0(i)*cl  dddT0(i)*c2  dhdTO(i), 

-1.5*dp(i)*cl  -1.5*dp(i)*c2  hp(i), 

1  0  0;0  1  0;0  0  l;-d(i)*c3  d(i)*c4  0]; 
A2=[-(rloc(i)+alt(i))*sm(lati(i))*cos(long(i)) 
-(rloc(i)+alt(i))*sin(lati(i))’'‘sin(long(i)) 

(rloc(i)+alt(i))*cos(lati(i)); 

-(rloc(i)+alt(i))*cos(lati(i))*sin(long(i)) 
(rloc(i)+alt(i))*cos(lati(i))*cos(long(i))  0; 
cos(lati(i))*cos(long(i))  cos(lati(i))''‘sin(long(i))  sin(lati(i))]; 
Atmp=Al*A2*r3(3*i-2;3*i,:);  %  A3=(delta  UEN  /  delta  xvz)=r3 

A=[A;(Atmp*dudUEN(i,:)’)';(Atmp*dvdUEN(i,:y)’];  %  total  A  matrix 

end 

AT==pinv(A);  Co  find  pseudoinverse  of  A 

dtp=AT*duv;  %  find  changes  to  tp 

tp=tp+dtp;  %  add  changes  to  tp 

ssel=sum(duv.'^2);  %  find  sum  of  the  squares  of  the  error 

dsse=abs(sse-ssel);  °o  difference  between  iterations 

if  ssel  <  sse 
sse=ssel; 
end 
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failindex=failindex+ 1 ; 
if  failindex>100 
dsse=0; 
end 
end 

%  Burnout  time  estimation 

Tlast=time(n); 

T  next=2  *time(n-2)-time(n-5), 
Tinax=TO+tmax; 
if  index  =1 
Tbo=T0+tmax; 
elseif  Tmax>=Tnext 
Tbo=Tlast-K).5*(Tnext-Tlast); 
else 

Tbo=Tlast+0. 5  *(T  max-Tlast); 
end 

tbo=Tbo-T0; 

%  Calculates  state  vector  at  burnout 


%  increment  counter 

%  if  statement  to  prevent  continuous  looping 
%  in  the  event  that  it  does  not  converge 


%  last  observation  time 
%  next  obser\  ation  if  TBM  was  still  burning 
%  max  obserxation  time  from  profile 

%  pick  bum  out  time  as  half  of  the  time 
%  between  Tlast  and  Tnext  or  Tlast  and 
®  b  Tniax  depending  on  relative 
"  o  magnitude  of  Tmax  and  Tnext 

“o  burnout  time-of-flight  (close  to  tmax) 


%  loop  to  next  tp  iteration 


d=(  1  - 1  5  *L)*(aO+al  •tbo+a2*tbo'^2+a3  *tbo^3+a4*tboM); 
ddot=(l-1.5*L)*(al+2*a2*tbo+3*a3*tbo^2+4*a4*tbo^3); 
h=(  1  +L)*(bO+b  1  *tbo+b2*tbo^2+b3  *tbo^3+b4*tboM); 
hdot=(  1  +L)*(b  1  +2*b2*tbo+3  *b3  *tbo^2+4*b4*tbo^3); 
theta=d/reff. 


®  o  distance 

horizontal  speed 
®  'o  height 
®/o  vertical  speed 
%  earth-central  anule 


latbo=pi/2-acos(cos(theta)*  sin(latO)+sin(theta)  *  cos(latO)*cos(hdgO));')  ogeod  t  i  c  1  at  i  t  u  d  e 
lonbo=lonO+asin(sin(theta)*sin(hdgO)/cos(latbo)); 
hbo=hO+h; 

Vbo=sqrt(ddot^2+hdot'^2), 
fpabo=atan(hdot/ddot); 


%  burnout  longitude 
%  burnout  altitude 
%  burnout  velocit\ 
%  flight  path  angle 


hdgbo=rem(asin(cos(latO)*sin(hdgO)/cos(latbo))+(2*pi),(2*pi));  °b  burnout  heading 


%  Calculates  ballistic  trajectorx  of  target  to  impact 

altbo=hO+hbo;  %  burnout  altitude 

gclbo=atan(((l-f)''2)*tan(latbo)),  ®  o  geocentric  latitude 

rlocbo=re*(l-f)/sqrt((l-f)^2*cos(gclbo)'^2+sm(gclbo)''2);  %  local  radius 

tgtposbo=[(rlocbo*cos(gclbo)+altbo*cos(latbo))*cos(lonbo) 
(rlocbo*cos(gclbo)+altbo*cos(latbo))*sin(lonbo) 
rlocbo*sin(gclbo)+altbo*sin(latbo)];  %  TBM  position  at  burnout 
rbo=norm(tgtposbo);  %  radius  at  burnout 

vo=rlocbo*we*cos(gclbo);  %  initial  inertial  velocity 

Vu=Vbo*sin(fpabo),  “  o  up  component  of  velocity 


99 


Ve=Vbo*cos(]^abo)*sin(hdgbo)+vo;  %  east  component  of  velocity 

Vn=Vbo*cos(^abo)*cos(hdgbo);  %  north  component  of  velocity’ 

Vin=sqit(Vu.'^2+Ve.^2+Vn.'^2);  %  inertial  speed 

fpain=asin(VuA^in);  %  inertial  flight  path  angle 

hdgin=Tem(pi/2-atan2(-Ve,Vin)+2*pi,2*pi);  %  inertial  heading 

Qin=Vin'^2*rbo/mu;  %  energy  parameter 

ein=sqrt(l+Qin*(Qin-2)’'‘cos(fpainy'2);  %  eccentricity  of  orbit 

PSI=2*acos((l-Qin*cos(fpain^2)/ein)+theta;  %ffeeflight  eca 

E=acos((ein-cos(PSI/2))/(l-dit*cos(PSI/2)));  %  eccentric  anomaly 

ain=rbo/(2-Qin);  %  semimajor  axis 

tfi=2*sqrt(ain^3/mu)*(pi-E+ein*sin(E));  %  time  of  free  flight 

gclim=asin(sin(gcIbo)*cos(PSI)+cos(gclbo)*sin(PSI)*cos(hdgin));  %geocentric  latitude 
latim=atan(tan(gclim)/(l-f)^2);  %  geodetic  latitude  of  impact 

flflon=acos((cos(PSI)-sin(gcIiin)*sin(gclbo))/(cos(gclim)*cos(gclbo)))-we*tfF; 
lomm=lonbo+fHon;  %  longitude  of  impact 

tini=2*tbo+tflF,  %  total  flight  time 

Tim=T0+tim,  %  time  of  day  of  impact 

data=[data;tp(l)  tp(2)  tp(4)  tp(3)  tp(5)  tp(6)  lonbo  latbo  hbo  fpabo  hdgbo  Vbo  Tim 
lonim  latim];  %  sa\  e  data  in  matrix 

end  %  end  of  numpoint  loop 


®  0  Plottinu  routines  and  data  anah  sis 


tllon=data(2,3  )*r2d; 
tllat=data(2,4)*r2d; 
llon=data(3  :index+l  ,3)*r2d; 
Ilat=data(3 :  index+ 1 ,4)’'‘r2d; 
mllon=mean(data(3  ;index+l  ,3))*r2d; 
inllat=mean(data(3  :index+ 1 ,4))*r2d; 
dllon=tllon-inllon; 
dUat=tlIat-mllat; 


°o  true  launch  lonuitude 
®  0  true  launch  latitude 
corrupt  launch  longitude 
%  coraipt  launch  latitude 
®'b  mean  of  corrupt  longitudes 
mean  of  corrupt  latitudes 
%  diff  between  mean  and  true  lion 
%  diff  betw  een  mean  and  true  Hat 


tblon=data(2,7)*r2d, 
tblat=data(2,8)*r2d; 
tbalt=data(2,9); 
blon=data(3 :  index+ 1 ,7)*r2d; 
blat=data(3  :index+ 1 ,8)*r2d; 
balt=data(3 ;  index+ 1 ,9); 
mblon=mean(data(3 :  index+1 ,7))*r2d; 
mblat=inean(data(3 :  index+ 1 ,8))*r2d; 
mbalt=mean(data(3 :  index+ 1 ,9)); 
dblon=tblon-mblon; 
dblat=tblat-mblat; 
dbalt=tbalt-mbalt; 


%  true  burnout  longitude 
%  true  burnout  latitude 
%  true  burnout  altitude 
%  corrupt  burnout  longitude 
%  corrupt  burnout  latitude 
“  o  corrupt  burnout  altitude 
%  mean  of  corrupt  bo  longitudes 
%  mean  of  corrupt  bo  latitudes 
%  mean  of  corrupt  bo  altitudes 
%  diffbetween  mean  and  true  blon 
%  diffbetween  mean  and  true  blat 
%  diffbetween  mean  and  true  blat 
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tilon=data(2, 14)*r2d; 
tilat=data(2, 1 5)*r2d, 

Uon=data(3 :  index+ 1 , 1 4)*r2d; 
ilat=data(3 :  index+ 1 , 1 5)*r2d, 
inilon=mean(data(3  :index+l,  14))*r2d; 
milat=niean(data(3 :  index+ 1 , 1 5))*r2d; 
dilon=tilon-milon; 
dilat=tilat-inilat; 

lpts=[llon-tllon  llat-tUat]; 
lk=cov(lpts); 

[lv,ld]=eig(lk); 

lab=2*sqrt(diag(ld)); 

la=lab(l); 

lb=lab(2); 

lth=pi/2-atan(lv(  1 , 1  )/lv(2, 1 )); 

bpts=[blon-tblon  blat-tblat  balt-tbalt]; 
[bv  1  ,bd  1  ]=eig(cov(bpts(: ,  1  ),bpts(:  ,2))), 
babc=2*sqrt(diag(bd  1 )); 
ba=babc(  1  ),bb=babc(2);bc=babc(3); 

ipts=[ilon-tilon  ilat-tilat]; 
ik=cov(ipts); 

[iv,id]=eig(ik), 

iab=2*sqrt(diag(id)); 

ia=iab(l), 

ib=iab(2); 

ith=pi/2-atan(iv(  1 , 1  )/iv(2, 1 )); 

%  Save  data  to  analyze 


%  taie  impact  longitude 
%  true  impact  latitude 
%  cornjpt  impact  longitude 
%  corrupt  impact  latitude 
%  mean  of  corrupt  im  longitudes 
%  mean  of  corrupt  im  latitudes 
%  diff  between  mean  and  true  ilon 
%  difT  bet  ween  mean  and  true  ilat 

%  center  data  about  launch  truth 
%  covariance  of  scatter 
%  eigen  value  &  vector  of  covariance 
%  dimensions  of  ellipse 
%  launch  ellipse  semimajor  axis 
%  launch  ellipse  semiminor  axis 
angle  to  rotate  ellipse 

°  0  center  data  about  burnout  truth 
®  0  eigen  val  and  vector  of  covariance 
°b  dimensions  of  ellipse 
"o  burnout  ellipse  semimajor  axis 

%  center  data  about  launch  truth 
®  o  co\  ariance  of  scatter 
%  eigen  value  &  vector  of  cox  arance 
°  o  dimensions  of  ellipse 
%  impact  ellipse  semimajor  axis 
%  impact  ellipse  semiminor  axis 
%  angle  to  rotate  ellipse 


data=[data;la  lb  1th  bal  bbl  bthl  ba2  bb2  bth2  ba3  bb3  bth3  ia  ib  ith;  mllon  mllat  mblon 
mblat  mbalt  milon  milat  dUon  dllat  dblon  dblat  dbalt  dilon  dilat  0]; 
ifii==0 

save  bOt  data  -ascii  -double 
elseifii=l 

save  b3t  data  -ascii  -double 
elseifii=2 

save  b6t  data  -ascii  -double 
elseifii==3 

save  b9t  data  -ascii  -double 
elseifii=4 

save  bl2t  data  -ascii  -double 
elseifii=5 
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save  bl  5t  data  -ascii  -double 
elseifii— 6 

save  bl8t  data  -ascii  -double 
elseifii=7 

save  b21t  data  -ascii  -double 
elseifii=8 

save  b24t  data  -ascii  -double 
elseif  ii=9 

save  bt27  data  -ascii  -double 
elseif  ii=  10 

save  bt30  data  -ascii  -double 
elseif  ii=  11 

save  bt33  data  -ascii  -double 
end 

clear  data 
end 
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